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ABSTRACT: The transport suspended sediment turbulent flow examined 
through the continuum theory mixtures. Turbulent averaging, applied the mixture 
momentum balances, yields correlation representing the effect the particles the 
drag associated with velocity fluctuations the fluid. This correlation balances the 
buoyant weight the sediment plane uniform mean flow. Thus, classical results 
for “turbulent particles may reinterpreted models for this 
correiation, and the source the diffusion process clearly identified view the 
dynamics fluid-particle interaction. This viewpoint leads multi-layer model for 
the mass transport turbulent boundary layer and results 
improvement the accuracy predicted concentration profiles steady plane flows. 


REFERENCE: McTigue, David F., “Mixture Theory for Suspended Sediment 
Transport,” Journal the Hydraulics Division, ASCE, Vol. 107, No. HY6, Proc. 
Paper 16313, June, 1981, pp. 659-673 
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KEY WORDS: Advection; Dispersion; Dyes; Hydrodynamics; Hydrologic 
models; Mississippi River; Mixing; Rivers; Simulation; Solutes; Transport 
theory; Water pollution; Water quality 


ABSTRACT: The transport (advection and dispersion) slug soluble 
suspendable material river analyzed two mathematical descriptions. The 
models are compared and are applied dye study 10.7-mile (17.2-km) long 
reach the Upper Mississippi River upstream from pool no.2. One model mixed 
cells-in-series (CIS) model and the other advection-dispersion (AD) model. Bulk 
parameters used both models are related each other and three sets field data. 
The CIS model simpler formulation and application than the model. Both 
models give very similar results for the field study. Both models show good agreement 
for nearly symmetrical time-concentration distribution, but give only fair description 
skewed time-concentration curves. explanation provided. 


REFERENCE: Stefan, Heinz G., and Demetracopoulos, Alexander C., 
Simulation Riverine Transport,” Journal the Hydraulics Division, ASCE, Vol. 
107, No. HY6, Proc. Paper 16349, June, 1981, pp. 675-697 
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KEY WORDS: Critical depth; Flow profiles; Fluid flow; Froude number; 
Gradually varied flow; Open channels; Specific energy; Subcritical flow 


ABSTRACT: compound-channel Froude number derived reflect correctly the 
occurrence minimum specific energy open channels with overbank flow. The 
limitations conventional Froude-number formulations are described and comparisons 
with the proposed compound-channel Froude number are presented. The results 
laboratory investigation into the occurrence two points minimum specific energy 
are presented, and this phenomenon can occur for certain compound-channel 
geometries predicted the compound-channel Froude number. The experimental 
investigation indicates that the upper point minimum specific energy can properly 
considered the limit subcritical flow. 


REFERENCE: Blalock, Merritt E., and Sturm, Terry W., Specific Energy 
Compound Open Channel,” Journal the Hydraulics Division, ASCE, Vol. 107, 
No. HY6, Proc. Paper 16292, June, 1981, pp. 699-717 


16299 DIFFUSION-WAVE FLOOD ROUTING CHANNEL 
NETWORKS 


KEY WORDS: Backwater; Channels (waterways); Flood routing; Floods; 
Hydrography; Networks; Open channels; Rivers; Unsteady flow; Water flow 


ABSTRACT: nonlinear diffusion-wave model developed for routing floods 
dendritic open-channel networks. The model considers backwater effects from both 
upstream and downstream ends channel. overlapping-segment technique 
adopted decompose the network. The nonlinear algebraic finite-difference equations 
are solved through special formulation tridiagonal banded matrices. Comparison 
the results the proposed model those the dynamic-wave and nonlinear 
kinematic wave models indicates that the proposed model nearly accurate as, and 
computationally much cheaper, than the dynamic-wave model. computationally 
inexpensive but more accurate than the nonlinear kinematic-wave model. 


REFERENCE: Akan, Ali Osman, and Yen, Ben Chie, “Diffusion-Wave Flood Routing 
Channel Networks,” Journal the Hydraulics Division, ASCE, Vol. 107, No. HY6, 
Proc. Paper 16299, June, 1981, pp. 719-732 


16307 MODELS SUSPENDED SEDIMENT DISPERSION 


KEY WORDS: Dispersion; Mathematical models; 
Probability theory; Sedimentation; Stochastic processes; Streams; Suspended 
sediments; Turbulent flow 


ABSTRACT: Existing probabilistic methods for predicting turbulent dispersion 
suspended sediments are reviewed, starting with the early theory particle diffusion 
continuous movements, and proceeding more recent random walk simulations. 
The evolutions these models are explored along with their limiting assumptions. The 
methods are basically equivalent. The random walk procedure not restricted 
principle stationary flows, but imposes probability distribution restrictions that are 
not present the continuous movements technique. basic problem common both 
techniques that detailed information needed about the turbulent features the 
carrying flow. The sediment diffusivity can obtained from either Lagrangian 
Eulerian information provided that the particle properties and turbulent velocity scales 
satisfy certain specific relationships. 


REFERENCE: Alonso, Carlos V., “Stochastic Models Suspended-Sediment 
Dispersion,” Journal the Hydraulics Division, ASCE, Vol. 107, No. HY6, Proc. 
Paper 16307, June, 1981, pp. 733-757 


16324 RUN HYDROGRAPHS FOR FLOODS 

KEY WORDS: Bivariate analysis; Damage; Dams; Floods; Hydrograph 

analysis; Hydrographs; Hydrology; Peak floods; Probability distribution 
functions; Spillways; Truncation errors; Volume 

ABSTRACT: bivariate lognormal probability distribution flood peaks and flood 
volumes describes families flood hydrographs, each family with desired return 
period. Truncated series floods are used estimate the underlying population 
parameters. The method provides more information floods than conventional 
methods, and the theory illustrated with worked example. Run hydrograph theory 
was developed and demonstrated for application fully-gaged catchments. The 
method allows for the synthesis the whole family design hydrographs for each 
desired return period. the development the theory, hydrographs were analyzed 
they occur the record; thus the sample hydrographs contains many composite 
hydrographs. This should result realistic design hydrographs for practical 
applications. the estimation the necessary parameters was demonstrated how 
straightforward can calculate the underlying probability distribution parameters 
from truncated sample. The run hydrograph technique may also amenable towards 
extrapolation ungaged catchments. 

REFERENCE: Hiemstra, Lourens A.V., and Francis, David M., “Run Hydrographs 
for Prediction Flood Journal the Hydraulics Division, ASCE, Vol. 
107, No. HY6, Proc. Paper 16324, June, 1981, pp. 759-775 
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MIXTURE THEORY FOR SUSPENDED 
SEDIMENT TRANSPORT 


David 


INTRODUCTION 


The transport fine particulate matter turbulent fluid flows fundamental 
importance broad range disciplines. Examples include the movement 
dust and pollutants atmospheric flows, sediment transport natural waters, 
and various industrial processes such the piping suspensions. 

steady mean rectilinear flow carrying particles denser than the fluid poses 
obvious problem: one must appeal some mechanism balance the buoyant 
weight the particles. was long ago suggested that this balancing mechanism 
could found diffusion’’ particles down mean concentration 
gradients. Indeed, results from diffusion models various forms compare quite 
well with both laboratory and field data for concentration profiles plane 
turbulent flows. 

The success such models has led considerable body work directed 
toward identifying the origin the process and toward generalizing 
the governing equations that more complex flows can analyzed. Most 
these efforts, however, assume priori gradient diffusion-type process, 
and yield little physical insight into the balance they attempt elucidate. 

This paper invokes the continuum theory mixtures examine the problem 
described the previous paragraphs. Mixture theory has been available for 
number years and has fact been applied the suspended sediment 
problem the past. However, the balance identified here has apparently been 
overlooked, and leads toa very simple and intuitively appealing result. Equilibrium 
struck between the buoyant weight the particles and correlation associated 
with the fluid drag due vertical velocity fluctuations. The classical 
term may interpreted model for this correlation. 

Having identified this correlation with the diffusion process, the scaling 
the diffusion coefficients clarified. two-layer model suggested, which 
the inner layer (near the bed) dominated the length scale the distance 
from the boundary, and the outer layer dominated the scale the flow 

Research Assoc., Dept. Earth and Planetary Sci., Massachusetts Inst. 
Tech., Cambridge, Mass. 02139. 

Note.—Discussion open until November 1981. extend the closing date one month, 
written request must filed with the Manager Technical and Professional Publications, 
ASCE. Manuscript was submitted for review for possible publication September 
1980. This paper part the Journal the Hydraulics Division, Proceedings the 
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depth. The predicted concentration profile, actually combination previously 
known results, agrees very well with available data. 


Previous Work 


Schmidt (18) first posed the governing equation for the suspended sediment 
distribution plane rectilinear flow noting that the net vertical flux 
particles must everywhere zero. The downward flux due gravitational 
settling must then exactly balanced upward diffusive flux, giving (see 


which concentration volume fraction solids; settling velocity 
particles; and diffusion coefficient for sediment. Eq. presented 


FIG. Sketch for Gravity Flow 


virtually all standard texts sedimentation (e.g., 16, 21), and various 
choices the coefficient, general taken function are explored. 
Considerable effort has been given developing more general theory that 
will yield Eq. special case. Several approaches are considered briefly 
below. 

Convection-Diffusion Theories.—Sediment dispersal turbulent flows having 
been observed exhibit diffusion-like characteristics, has often been suggested 
that classical convection-diffusion model for the concentration appropriate 
generalization. This can motivated the following manner. The mass balance 
for the fluid, constant density, p’, given 


Fig. 1): 
(1) 
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which fluid velocity. Introducing the Reynolds decompositions, 
indicate fluctuating quantities, and averaging gives 


t 


Boussinesq model for the correlation now introduced, such that 


which second-order tensor diffusion coefficients, Eq. becomes 
the classical convection-diffusion equation: 

Eq. does not reduce Eq. for plane rectilinear flow, although that assertion 
has erroneously been made (9). Furthermore, Eq. motivated entirely 
kinematic arguments, and provides insight into the dynamic process sediment 
dispersal. 

Two-Phase Flow and co-workers (20) have developed 
extensive body theory and analysis based part the dynamics the 
particulate phase. This theory able reproduce the balance given Eq. 
but, its present form, the theory exhibits several shortcomings. 

First, the gradient diffusion term introduced by, for example, Soo and Tung 
(20) again has dynamic significance. arises through the Reynolds decomposi- 
tion and averaging applied the particulate phase mass balance (19), 
fashion similar that shown Eqs. further noted here that, although 
the averaging scheme applied the mass balance obtain the diffusion 
term, not then applied the momentum balance. Thus, other correlations 
fluctuating quantities that should accounted for are neglected. Also, this 
leads situation where averaged velocities are equated with exact velocities. 
the averaging carried out consistently throughout the derivation the 
two-phase theory, results very similar those that follow can obtained. 

Finally, the two-phase theory leads second-order equation for the 
concentration plane rectilinear flow, the first integral which gives Eq. 
when the integration constant taken zero. This introduces further 
ambiguity the problem already present determining appropriate condition 
for even the first-order equation given Eq. 


Mixture 


The approach adopted here grounded the formalism the continuum 
theory mixtures [see Atkin and Craine (1,2) for general review]. few 
studies particulate transport turbulent flows have been made view 
mixture theory, but have led significantly different results Comparisons 
will noted where appropriate the following developments. 

The theory assumes that the fluid and the dispersed solid particles can each 
treated continuous media. For fluid and solid densities respectively, 


7 
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and particle volume fraction (concentration) the bulk densities the fluid 
and particulate phases the mixture are and respectively. The 
mass balance equations, for constant and are then 


which velocity dispersed solid phase. The momentum balances are 


(11): 


a 
t 


t 


which fluid phase stress; dispersed solid phase stress; acceleration 
gravity; and force the fluid the particulate phase per unit volume. 
Note that the interaction term, vanishes when and are added 
obtain momentum balance for the mixture. Also, the limit vanishes, 
and reduce properly the momentum balance for single fluid. 
Closure the system Eqs. 6-9 obtained posing constitutive equations 
for the stresses, and and the interaction force, For simplicity, the 
independent constitutive variables assumed here are 


both which are properly invariant (14). Neglect any dependence the 
rate deformation the fluid, the rate deformation the particulate 
phase, and the relative rotation rate anticipates that concern limited here 
flows which inertia dominates the flow behavior. Thus, the detailed form 
the rate-dependent terms the constitutive equations unimportant, 
assumption that clearly valid the limit very small concentrations 
and high mean flow Reynolds numbers. The neglected dissipative stress 
terms are crucial, course, energy balance considerations (5), which are 
not addressed here. 


The linear, isotropic constitutive equations based the assumptions Eq. 


which and the fluid and dispersed solid phase pressures, respectively; 
and and general, scalar functions the concentration, Thus, each 
phase modeled ideal fluid. 


assumed that the pressure each phase identical; that is, 


are 
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(7,11). Thus, simply the mean normal stress the mixture, (1/3) 
are readily identifiable with well-known fluid-particle interaction forces. The 
term the relative velocity, the drag force. its linear form, 
the drag term limited the Stokes range (small particle Reynolds numbers), 
and, for noninteracting spherical particles, given 


which particle radius; and fluid viscosity. The validity the drag 
term may easily extended larger particle Reynolds numbers allowing 
nonlinear terms Eq. 13, which case can function the invariant 
(20). The second term Eq. identified with the force the 
fluid pressure the dispersed solid phase, giving (6) 


The form the term suggestive gradient diffusion process, 
and indeed, diffusivity, might included represent Brownian 
diffusion. Peddieson (15) has suggested that such term may also represent 
diffusion’’ the particles. However, simply inserting turbulent 
diffusivity into the constitutive model somewhat hoc, insofar such 
process cannot readily identified the exact equations motion. Reynolds 
decomposition and averaging the Brownian diffusion term does not result 
additional ‘‘turbulent term, and the effect Brownian motion 
large particles exceedingly small. 

Collecting the results Eqs. 8-15 


a? 


au’ 


Ss Op 


16, and give complete set eight equations for the eight unknowns, 
u’, and Identical equations motion have been used Drew (4) 
and derived Hinze (10), neglecting viscous terms and virtual mass effects. 


AND AVERAGING 


The decomposition and averaging scheme adopted here from Drew (4), 
and only briefly outlined the following. The velocities and concentration 
momenta, such that 


Thus, the average the velocity fluctuations does not vanish. Rather, 
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cu’’ and cu” =0. originally developed for treating compressible 
fluid flows, yields only one correlation from the inertial terms Eqs. 
17. Applying these definitions Eqs. 16, and gives the mass and 
momentum balances for the fluid and particulate phases: 


a 


were obtained Drew (4). 


Consider now the classical problem steady plane gravity flow with free 
surface (Fig. 1). Solutions are sought for the flow when 
and ¢(y). Eqs. and are then identically satisfied, and and 


0 
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Adding the two x-momentum equations, Eqs. and 26: 


When and are the same order, the case for, say, quartz particles 
water, and assuming that the fluid and particle velocity fluctuations are 
the same order, Eq. may approximated for 


the classical balance between the gravitational driving force and the gradient 
the Reynolds stress pure fluid. Note that Eq. suggests that the 
characteristic scale the velocity fluctuations the order (gh sin 
which the flow depth, leading the usual definition shear velocity; 
Eq. 29, course, requires the adoption model for the correlation 
The most common choice mixing-length model, which leads the present 
case the familiar logarithmic profile for (16), and requires 
further consideration here. 

Now consider the y-momentum balances, Eqs. and 27. Eliminating the 
mean pressure gradient between the two equations and again approximating 
for small concentrations and small slopes (cos gives 

7 9 


+ p’—(cv” — fp’ —(v”’ ) 


The second term the buoyant weight the dispersed particulate phase, and 
the term terms that balance this are sought. order examine this balance, 
Eq. nondimensionalized assuming that scales with characteristic 
concentration c,, v’’ and with (see aforementioned identity), and 
with 1/h. Rearranging that the buoyant weight term order unity, Eq. 
becomes 


characteristic mean flow velocity; and all variables 
are now dimensionless and order unity. Note that the group 
the square the Froude number. 
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Consider typical wide stream carrying very fine quartz sand 5.0 
bottom. For slope and depth 1.0 (39 in.), the mean flow 
velocity about 1.0 m/s (3.3 fps) and found that the second term 
Eq. order unity, while the terms brackets are order 10~*. Thus, 
returning dimensional variables, clear that the buoyant weight the 


particles Eq. balanced the drag associated with the vertical velocity 
fluctuations the fluid: 


Note that Drew (4) examined the identical problem view the same governing 
equations, but discarded the term early his analysis. was thus 
forced find the balance with the Reynolds stress terms appearing Eq. 
30. 

The foregoing arguments assume that the characteristic length scale the 
flow the total depth. Near the bed, however, the relevant scale 
the distance from the bed. This suggests, based 29, that the velocity 
fluctuations scale with very near the bed. The dimensionless group 
scaling the second term Eq. then becomes Thus, the 
terms brackets Eq. may dominate over the second term region 
very close the lower boundary. For the sample parameters cited previously, 
however, these additional terms would important only vanishingly thin 
layer, and the balance given Eq. appears valid over virtually the 
entire flow depth. 


Dividing Eq. recovers the balance terms the flux 
settling particles, 


and comparison with Eq. shows immediately that the classical 
may regarded model for the correlation .Similar balances 
have been noted before, based simply the kinematics particle fluxes. 
For example, both the ASCE (21, pp. 69-71) and Lumley (12, 309) obtain 
expressions like Eq. 33, but with v’’ replaced nonetheless particularly 
satisfying discover this balance the dynamics the fluid-particle interactions. 

The interpretation Eq. straightforward. Due viscous drag the 
particles, vertical turbulent motions fluid parcels carry sediment. This process 
tends smooth out gradients, thus appearing 
The correlation Eq. clearly nonzero the presence 
vertical gradient the mean concentration. For increasing downward, 
upward (+) fluid velocity fluctuation carries with concentration larger 


(+) than the local mean, and downward (—) velocity fluctuation carries 
concentration deficit (—). 


view the foregoing results, various expressions for 
sediment are interpreted here models for the correlation This 
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term arises from Reynolds decomposition and averaging applied the fluid-par- 
ticle interaction force appearing the mixture momentum balances (Eqs. 
and 9). Although originally motivated different ways, number well-known 
theories for diffusion’’ exist, and three are reviewed briefly here 
light the preceding reinterpretation. All are variations the Boussinesq 
model for the correlation given Eq. noting that the averaging scheme 
now that defined Eqs. and 19. 
The simplest expression for the diffusion coefficient tensor 


which scalar. Then Eq. becomes 


and substitution into Eq. for the plane uniform flow retrieves the classical 
balance Eq. inspection Eq. 35, can anticipated that the diffusivity, 
€,, scales with the product the characteristic velocity and length scales 
the flow. direct analogy with models for momentum transport turbulent 
boundary layer, expected that the concentration profile will exhibit 
inner and outer region because the dominance different length scales 
near the bed and far away. 

The outer region commonly taken (8) include about 80% the boundary 
layer (in the present case, the boundary layer occupies the entire flow depth). 
this region, the characteristic velocity scale and the dominant length 
scale the flow depth. Thus 


which dimensionless constant order unity. For the outer region, 
then, constant, and Eq. easily solved, yielding the familiar result: 


The boundary condition has been used here, where some arbitrary 
height above the bed. Classically, this arbitrariness due lack anything 
the theory that would suggest the appropriate condition true boundary. 
the context the multilayer model presented here, the problem identifying 
appropriate condition the near-bed concentration resides the inner layer, 
while Eq. obtained matching with the inner solution some 
level, say 0.2h. 

Eq. shown Fig. fitted data from Montes and Ippen (13) for 
1.4 m/s (4.6 fps), 5.7 (2.2 Note that the particle Reynolds 
number falls slightly outside the Stokes range, that the drag term linear 
the theory not strictly valid. Nonetheless, the agreement 
excellent the outer 80% the flow, while, anticipated, Eq. fails near 
the bed. value 0.11, which the correct order magnitude, 
obtained for this example. 


The constant given Eq. implies isotropic, homogeneous turbulence, 
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Equation 
K,20.11 


run 37, 
Montes 
Ippen 


FIG. 2.—Concentration Profile for Outer Region 


Equation 


run 37, 
Montes 


FIG. 3.—Concentration Profile for Inner Region 
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and must clearly break down near solid boundaries where the correlation 
must vanish vanishes. Just the usual arguments for momentum 
transport, the appropriate length scale near the bed while the velocity 
scale u,, giving for the inner layer 


which again dimensionless constant order unity. This model has 
been suggested Batchelor (3) and Lumley (12), and gives, with Eq. 


which the condition has again been imposed. This result appealing 
that the concentration becomes more uniform for increasing bottom shear 
Stress, and for decreasing settling velocity. Note that Eq. gives 
Eq. shown Fig. for comparison with the same experimental 
data Montes and Ippen (13), redrawn for ease interpretation log-log 
plot. value 0.35 obtained for the fit shown. The agreement 
again excellent, this case for the inner 50% the flow. Thus, dividing 


Equation 
0.44 


run 37, 
Montes 
Ippen 


FIG. 4.—Concentration Profile from Rouse (17) Theory 


h-y 
2 
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the flow into inner and outer layers characterized different length scales, 
very good predictions are obtained for each region, and they can matched 
their broad overlap region. 

Finally, third expression for was given Rouse (17), who supposed 
that the Eq. must similar the eddy diffusivity 
introduced considering the fluid velocity field. Thus, taking 


the resulting concentration profile, from Eq. 33, 


h -y 
(41) 


Eq. shown Fig. again with the same data for comparison, the fit 
shown giving 0.44. The model Eq. cannot motivated unambiguously 
from dimensional arguments Eqs. and 38. However, recognizing that 
both and are important length scales, Eq. might taken the simplest 
model that ensures that the correlation vanishes both the bed and the 
free surface. 

The widely accepted result Eq. enjoys the advantage being 
solution, valid over the entire flow depth. However, from what known about 
momentum transport turbulent boundary layers, there every reason 
expect that mass transport, too, can more accurately represented multilayer 


model. Available data, including those shown Figs. 2-4, indeed bear this 
out. 


Summary 


The transport sediment turbulent flow over flat bottom has been 
examined view mixture theory, which models both the fluid and particulate 
phases continua. The momentum balances for the two phases are coupled 
through interaction terms that include pressure gradient and drag forces. 
decomposition and averaging scheme accounts for turbulent fluctuations the 
velocities and concentrations. The resulting equations motion have been 
presented previously Drew (4). 

For plane uniform flow down slope low concentrations, the sediment 
does not affect the bulk flow behavior, and the usual balance between the 
gravitational driving force and the gradient the fluid Reynolds stress obtained. 
The cross-stream momentum balances yield equation that includes the buoyant 
weight the particles, and scaling arguments suggest balance with the 
plete the identification this correlation with the process 
and classical results are retrieved. 

The concept balance between the particle flux due gravitational settling 
and diffusive flux due turbulence has long been established. number 
studies have proposed generalizations this concept, but many simply assume 
priori the existence Fickian diffusion-like process. Such theories can 
only model observation, rather than offering any insight into its source. 
Other studies identify correlation similar that Eq. from purely kinematic 
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arguments, again providing information the dynamics the dispersal 
process. 

The theory presented here has identified the process terms 
correlation that arises from the decomposition and averaging the 
fluid-particle interaction terms the momentum balances. The result simple 
and intuitively appealing. The correlation represents the effect the drag 
the particles associated with turbulent velocity fluctuations the fluid. the 
case plane flow over flat bed the presence mean concentration 
gradient, clear that the concentration and vertical velocity fluctuations 
are highly correlated. The result will appear diffusion-like process, 
tending obliterate the mean concentration gradient. 

Thus, mixture theory offers clear dynamic interpretation the 
sediment. The diffusion process, successfully treated Rouse 
(17) and others, may then viewed model for the correlation identified 
here. 

Review previous representations (Eqs. 36, 38, 40) the 
diffusivity,’’ €,, light the present theory suggests that more accurate results 
can obtained combining them multilayer model. Just momentum 
transport processes turbulent boundary layer, mass transport can divided 
into inner and outer regions. The outer part the flow dominated the 
length scale while the inner part, near the bed, dominated Joining 
the solutions for these two layers overlap region yields results that compare 
more favorably available data than the widely accepted Rouse equation (Eq. 

For clarity and simplicity, the data presented here are confined single 
example from Montes and Ippen (13). However, examination other data 
consistently reinforces the success the layered model considered here. 
This has been done, e.g., for the data presented the ASCE (21, 80) 
support the Rouse equation (Eq. 41). The prediction for the concentration 
the outer layer (Eq. 37) significantly larger than few measured values 
very near free surface. This might anticipated, since the correlation 
does not vanish the free surface under the assumption Eq. 36, while 
Rouse’s formulation (Eq. 40) assures this. the context the multilayer model 
presented here, this difficulty can easily accommodated introducing 
third layer dominated the length scale 

addition recovering well-known classical results, the present theory is, 
from the outset, fully three-dimensional (Eqs. 20-23). further advantage 
that, the level Eqs. and 17, the structure the theory completely 
general. Thus, directions for future refinement are immediately apparent. First, 
the constitutive equations used are based the greatest possible simplifications. 
Considerable work needed develop appropriate forms for the stresses and 
interaction forces that include the rate-dependent constitutive variables and that 
are valid higher concentrations. The second major area for further research 
clearly developing appropriate scaling relations and models for the correla- 
tions appearing the theory. 
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concentration (volume fraction solids); 
reference concentration; 

Brownian diffusivity Boltzmann’s constant, 
absolute temperature); 

force fluid particulate phase per unit volume; 
scalar functions 

acceleration gravity; 

flow depth; 

pressure; 

fluid phase pressure; 

dispersed particulate phase pressure; 
time; 

fluid phase stress; 

dispersed particulate phase stress; 

fluid phase velocity; 

characteristic mean flow velocity; 
single-particle settling velocity; 

dispersed particulate phase velocity; 
shear velocity; 

components 

Cartesian coordinates; 

bed slope; 

diffusion coefficient tensor; 


sediment diffusivity; 

dimensionless constants order unity; 
arbitrary distance above bed; 

fluid viscosity; 

fluid density; 

solid density; and 

bottom shear stress. 


Symbols 
overbar averaged quantity; and 
prime fluctuating quantity. 
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CELLS-IN-SERIES SIMULATION 
RIVERINE TRANSPORT 


INTRODUCTION 


cells-in-series (CIS) formulation simulate riverine transport dissolved 
material water quality models will presented. The characteristics the 
model will compared one using the conventional one-dimensional advection- 
dispersion (AD) equation, and the advantages and disadvantages the CIS 
model will illustrated. Both models will applied 10.7-mile (17.2-km) 
reach the Upper Mississippi River which dye tracing experiments have 
been conducted the United States Geological Survey (USGS). 

CIS model for the transport conservative material river was 
previously presented Banks basic description CIS modeling for 
chemical reactors can found Levenspiel (14). 

Major discrepancies between observed dispersion rivers and predictions 
with the classical equation have been acknowledged and explained the 
literature (16,20). Some investigators have therefore proposed other models 
predict mass transport natural streams and rivers. Examples are those given 
Thackston and Schnelle (24), Jobson and Yotsukura (12), McQuivey and 
Keefer (17,18), Pedersen (19), Sabol and Nordin (20), Beltaos (2), and Liu (16). 
Some the alternative models are expansions the equation and include 
entrapment dead zones (19,24,27), multidimensionality (12), time-depen- 
dent (or distance) dispersion coefficients (16). Models which not use the 
equation all emphasize convective transport (18) combination 
convective transport and stochastic processes (20). The convective model (18) 
was shown give very good results for the first one-sixth the convective 
period defined Fischer (5). 

The approach taken this paper more radical deviation from the classical 
equation. The CIS formulation used herein has been recognized 
alternative way for mass transport modeling essentially one-dimensional 
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systems (14). Only the fundamental aspects the CIS theory and its relation 
models will considered herein. The question how relate the 
CIS model parameter, (number cells), the hydraulic characteristics 
stream river will not covered detail. 

The physical process mass transport and dispersion natural river 
quite complex. Any mathematical description natural systems requires simpli- 
fications and abstractions. The search for simple yet appropriate description 
has led the CIS model which has basic equation that ordinary first-order 
differential equation instead partial differential equation, used 
models. This advantageous when the basic transport equation used 
water quality model for extensive river reach, branching system, 
the case shallow river impoundments, estuarine networks, ship canals, 
etc. 

has also been suggested reviewer that CIS simulation may applied 
approximation river reaches where the one-dimensional equation 
cannot used because the initial convective period (5) has not been exceeded. 
that case, two-dimensional model should really used. However, the 
complexity the river geometry the time and expense not justify 
development two-dimensional model for such river reach, has been 


suggested that CIS model with site-specific parameters may provide approximate 
answers. 


Transport River 


The following presentation will consider the transport slug material 
river because illustrates the riverine transport features which are 
interest water quality modeling. idealized but typical concentration 
versus time record which describes the passage slug dissolved fine 
suspended material two consecutive stations river shown Fig. 
Without specifying the mass concentration functions, C(t), analytically, some 
important bulk features these functions can described the moments 
C(t). shall consider the first four moments C(t) because they suffice 
for fairly accurate summary the C(t) curves for water quality modeling 
purposes. 

The four moments interest are: 


Cdt (zero moment) 


Ctdt (first moment) 


+o 


(second moment) 


o 


(third moment) 


Each the four moments related descriptive particular feature 
characteristic riverine transport. 


au 
v=o (3) 
(4) 
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Parameter M/Q, which the steady-state volumetric river flow 
rate; and the total mass material the slug. applied two stations, 
and respectively, along stream, (A, A,) represents the loss 
material between the two stations. Such loss may produced chemical 
biological reactions, settling, evaporation, etc. the mass conservative, 
concentration versus time function should give the same value all 
stations along river reach. 

The first moment, specifies the location (in time) the centroid the 
slug mass. differential, between two stations and the 
mean advective travel time residence time the material the river reach 
between these two stations. 

The variance, reflects the longitudinal spread material along the river. 
measure the width (in time) the two curves, C(t), shown Fig. 
but does not account for the lack symmetry about the mean travel time, 


STATION 


STATION 


TIME 


FIG. 1.—Schematic Concentration Versus Time Curves for Slug Dissolved (or 
Suspended) Material Traveling Past Two River Sampling Stations 


longitudinal dispersion, which turn depends such flow characteristics 
velocity distributions and turbulence. 

Finally, the third moment skewness, describes the asymmetry the 
curve, C(t). This asymmetry can caused the velocity distribution 
real river entrapment materials, separated flow regions and backwaters 
and subsequent release. Entrapment usually causes large degree skewness. 
The last three the aforementioned features are functions river morphology 
and flow rate, while the first (zero moment) simply the mass, divided 
the discharge, 

can shown how CIS model, and comparison, the equation, 
describe the aforementioned four features. 


CIS 


CIS model assumes that river composed sequence elements 
equal volume, called cells. Complete mixing occurs each cell. The general 
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KVC 


time; W(t) external mass input rate into the cell, including that advected 
into the cell from the neighboring cell upstream, and that through direct injection; 
the volumetric flow rate out the cell; and kinetic rate material 


For impulse input, characteristic for slug release, W(t) 
which the Dirac Delta function, such that when 
and 


The concentration the first cell 


and the concentration for the kth cell 


y* 


C,(t) 


The four transport characteristics defined Eqs. can now computed. 
the process, the gamma function 


integer, will used. The material conservation coefficient found 


The mean and the variance the gamma-type concentration distribution are 


1)! 
a 
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13-16 are simple functions the number cells used. This characteristic 
the CIS model more directly shown Fig. which the dimensionless 
concentration, the effluent from river reach constant volume, 
V,, has been plotted versus dimensionless time, Note that kV, 
which the number cells the reach; and the volume cell. 
The normalized form Eq. 


> 
2.0 


NORMALIZED TIME 


FIG. (Normalized) Concentration Versus Time Curves for Different 
Numbers Cells River Reach (Eq. 17) 
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Ki t k* 


Fig. was plotted for 

The loss mass, the mean residence time, and the changes variance 
and skewness between two stations and can expressed simple functions 


Because the CIS model simple model, several the relationships are linear 
with the number cells k,_,. 


The one-dimensional transport model for stream represented the 
relationship 


which the mean flow velocity the longitudinal dispersion 
coefficient; and the distance along the stream. The longitudinal dispersion 
coefficient, introduced analogy with Fick’s law diffusion, although 
the riverine dispersion process much more complex than, and has little physical 
resemblance to, diffusion process. extensive amount literature exists 
the relationship between the dispersion coefficient and stream and hydrologic 
stream characteristics. Some generally pertinent but not exclusive references 
are Refs. 15, 23, and 28. More recently, the effect temporary entrapment 
materials separated flow regions river systems has been studied (19, 
20, 25, and 27). 


The solution Eq. for impulse input tracer described earlier 
and shown Fig. 


4Dt 
the river flow rate, and the total volume the river reach downstream 


from the injection point, V,, are introduced, the equation can rewritten 
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This equation gives the concentration tracer function dimensionless 
time. The time normalized the mean residence time, The parameter 
the inverse Peclet Number, Pe. The skewness the C(t) curve 
increases with D/Ux and Kt. The parameter can also expressed 
for practical evaluations. For and very small values 
(large Peclet numbers), the concentration distribution approaches Gaussian 
function. The skewness predicted Eq. usually much less than that 
observed river and stream and expanded versions Eq. with additional 
Eqs. 16, 19, and 24) are required account for the real skewness. 

For conservative tracer 0), Eq. can used evaluate the transport 
parameters introduced 1-4. Eq. gives 


Finally, Eq. gives 


The last three parameters are related the Peclet number the flow. 

25-28 are also obtained provided that Eq. 24, multiplied 
constant time, e.g., rather than instantaneous time, The foregoing 
analysis therefore valid for all dissolved materials the water (conservative 
well nonconservative), provided that the kinetic decay rates are low. 


The mean travel times, variances, and skewnesses the CIS model expressed 
Eq. 14-16, respectively, can compared those the model, expressed 
Eqs. 26-28, respectively. relate the equations one needs replace the 
total river reach volume addition, (Q/V) and for 
conservative tracer V,/Q k/a. comparison the two sets, 
Eqs. 14-16 and Eqs. 26-28 shows that the expressions for and 
derived for the model degenerate into those the CIS model the tracer 
conservative and the Peclet number sufficiently large. Otherwise, 
the two models provide somewhat different results. 
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C(t) ... (24) 
According Levenspiel and Smith (13), Eqs. and give 
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given river reach volume, V,, and given river flow, the 
mean travel time independent the subdivision (number cells) the 
CIS model (Eq. 14), because k/a V,/Q. slightly dependent the 
dispersion coefficient, the models (Eq. 26). 

The spread the arrival times the dissolved (or suspended) material 
any one station related the number cells into which the river reach 
subdivided the CIS model. larger number cells provides less dispersion, 
shown Fig. The longitudinal dispersion coefficient, the model 
and the number cells, the CIS model therefore describe similar transport 
features and are related each other through 


PIGS EYE 
LAKE 


STATION 


STATION 


STATION 


STATION 


STATION 
STATION 


FIG. 3.—Upper Mississippi River Reach Downstream from St. Paul 


/ 
STATION 
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The number cells also affects the skewness the C(t) curve. does 
the longitudinal dispersion coefficient the models. Typically, the CIS 
model can provide larger degree skewness than the model, which 
often more realistic. both the CIS and the models, skewness not 
described independent parameter but determined the choice 

Other parameters can also selected describe the riverine transport features 
particular river reach, e.g., the maximum concentration the C(t) curve. 
For the CIS model found 


and for the model 


with the times arrival the maximum concentrations given 


Qa 


and 


for the CIS and the model, respectively. Eqs. and are for 


July 1976, May 1977, and October 1977, the USGS St. Paul District traced 
the passage dye releases the Upper Mississippi the vicinity Minneapo- 


TABLE 1.—Location Sampling Sites Upper Mississippi River 


Length 
USGS site Station Miles above subreach, 
number number Location Ohio River kilometers 


(1) (2) (3) (4) 


Pigs Eye Metropolitan Waste 
Treatment Plant 

1-494 Bridge 

County Road Bridge 

Shieley Commercial Dock 

Spring Lake Inlet Channel (July 
1976) 

Main channel Lower Gray 
Cloud Island (May 1977) 

Downstream from Spring Lake 
Entrance (October 1977) 


k-1 
a 
(33) 
(5) 
7.6 835.2 
832.5 
830.3 3.5 
826.5 6.1 
823.8 4.3 
823.8 
823.0 1.3 


JUNE 1981 HY6 
TABLE 2.—Upper Mississippi Transport Parameters 


July 1976° 
Station 
Parameter 


(1) 
kilograms 
kilograms 
hours 
square hours 
S(-) 
hours 


3,260 cfs (92.32 s); 52.19 kg. 


TABLE 3.—Upper Mississippi Transport Parameters Calculated 


July 1976 
Station 
Parameter 


(1) 


k(-) 

integer 

per hour 

per hour 

K/a 

cubic meters) 

milligrams per liter 

hours 

kV, 10° cubic feet 10° 
cubic meters) 


“The average values for the three surveys are 664 10° (183 10° 
Note: Values M,, and are the same Table 


Paul. The procedures and basic data were described Ref. 26. Rhodamine 
was injected three different locations the river upstream from Station 
(Fig. and traced downstream. The river reach from mile 835.2 mile 
823.8 downstream from the Metropolitan Waste Treatment Plant was covered 
all three experiments. (The mileage with reference the mouth the 
Ohio River.) For orientation, the Wabasha Bridge near downtown St. Paul 
river mile 839.5; Lock and Dam No. river mile 815.4. Fig. shows 
the Mississippi River reach interest this study and the locations the 
sampling stations. Table provides some additional information the sampling 
stations. 


(5) 
36.66 36.01 39.50 39.35 
0.702 0.690 0.757 0.754 
27.2 41.6 62.9 84.1 
75.69 124.2 213.3 226.5 
0.616 1.009 1.147 0.651 
6.0 4.4 3.7 3.45 
18.0 32.0 60.0 82.0 
9.75 13.94 18.55 31.19 
10 14 19 31 
0.359 0.335 0.295 0.371 
0.013 0.009 0.005 0.003 
0.965 0.972 0.984 0.991 
33.9 36.0 40.5 31.9 
(9.6) (10.2) (11.5) (9.03) 
0.641 0.536 0.464 0.358 
6.7 4.2 3.1 
24.4 38.6 59.5 81.4 
331 502 751 995 
(93.7) (142) (213) (282) 
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from Dye Release Data for Three Field Surveys 


May 1977° October 1977° 
Station Station 


from CIS Model Section for Three Field Surveys 


May 1977 October 1977 
Station Station 


600 10° (136 10° and 654 10° (187 10° respectively. 


Riverine transport parameters presented the preceding sections were calcu- 
lated from the dye concentration versus time data measured each station. 
The numerical results are presented Table The CIS model and the 
model were applied the field data follows. 


CIS Model 

Analysis Data River Cross Section.—If only river flow and dye study 
data but river morphologic data are used, the observed values 
and Table can used find the time constant, the number 
cells, and the kinetic rate, section section. For this purpose 
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(6) (7) (11) 
45.95 43.68 45.04 41.72 27.24 25.43 
0.921 0.875 0.902 0.836 0.632 0.590 
11.9 20.2 34.4 47.4 9.58 28.3 
14.4 30.1 37.3 66.4 3.64 34.3 
2.46 1.94 1.38 1.31 1.61 1.48 
22.0 14.4 9.4 6.4 10.5 4.0 
10.3 17.0 30.0 42.0 8.0 24.8 
(6) (7) (11) 
9.83 13.62 31.69 33.84 25.3 23.31 
0.825 0.673 0.922 0.714 2.64 0.824 
0.007 0.007 0.003 0.004 0.049 0.019 
0.991 0.990 0.997 0.995 0.981 0.977 
20.6 25.2 18.3 23.7 13.2 42.4 
(5.8) (7.14) (5.18) (6.71) (3.74) (12.0) 
0.638 0.542 0.355 0.344 0.398 0.414 
10.8 8.1 6.5 7.9 1.9 
10.7 18.8 33.3 46.0 9.2 27.1 
202 343 580 802 334 988 
(57) (97) (164) (227) (95) (280) 
2 
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13, 14, and are rewritten the form: 


The results such computations for the Upper Mississippi River are shown 
Table Concentrations down 0.1 were considered the receding 
limb the measured curves. With these parameters, other values such 
and can calculated using Eqs. 30, and 32, respectively. 

The kinetic rates, computed from Eq. were found virtually identical 


with those obtained the conventional method which measured values 
and are fitted equation 


which the recovery ratio the tracer. 
Analysis Data River Stations and the Mississippi River 


are chosen, e.g., the pertinent relationship derived from 13, 14, and 
are 


Measured 
Model 


Concentration 


12 13 16 17 18 
Time (days) 

July 1976 


FIG. 4.—Measured and Predicted Model) Concentration Versus Time and Station 
1-Station Dye Release Study, July 1976 
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Measured 
Model 


Concentration 


Time (days) 
May 


FIG. 5.—Measured and Predicted (CIS Model) Concentration Versus Time Station 
1-Station Dye Release Study, May 1977 


Meosured 
CiS 


Concentration 


NWA UD N O 


Time (days) 
Oct. 1977 


FIG. 6.—Measured and Predicted (CIS Model) Concentration Versus Time Station 
1-Station Dye Release Study, October 1977 
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TABLE 4.—Upper Mississippi Transport Parameters Calculated 


July 1976 
Station 


Parameter 


(1) (2) (3) 
V,/Q, hours 24.7 38.8 59.7 81.5 
per hour 0.00205 0.00093 0.00045 0.000196 
feet (meters) 14,250 25,870 45,930 60,190 
(4,340) (7,880) (14,000) (18,340) 
square feet per second 115.7 172.9 263.7 197.2 
(square meters per second) (10.75) (16.06) (24.50) (18.33) 
10° cubic feet 289.9 455.4 700.6 956.5 
10° cubic meters) (82.0) (128.9) (198.4) (270.8) 
per liter 5.67 4.17 3.33 2.84 
hours 37.4 58.1 80.2 
per foot (per meter) 0.057 0.057 0.057 0.057 
(0.189) (0.189) (0.189) (0.189) 


Note: Values and are the same Table 


M 
4 


The variable the number mixed cells between Stations and 
k,. Numerical values are not shown for lack space. The results 
were comparable those obtained analysis cross section Table 
Predicted Curves.—With the parameters and established, one 
may proceed calculate C(t) curves and compare them with the actually 
observed ones. This was done for all three field survey dates and each station. 
Predicted (CIS model) and measured concentration versus time distributions 
station are shown Figs. and The July 1976 concentrations 
Fig. and the October 1977 concentrations Fig. are reasonably 
well-predicted, except for some irregular features Fig. which can neither 
explained nor modeled. The predictions Fig. are good because the parameters 
and were selected station station from Table The May 1977 
data Fig. produce apparent poorer match between predicted and measured 
concentration versus time curves. particular, the maximum concentration 
poorly-predicted. The same discrepancy will seen occur for the 
model (see later figure) and explanation will provided later. 


Model 
Analysis Data River Cross Section.—From Eqs. and 27, relationships 
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from Model Section for Three Field Surveys 


May 1977 October 1977 
Station Station 


for the dispersion coefficient, and the total hydraulic residence time, 


functions the measured parameters, and can developed. 
found that 


Numerical results are shown Table 


Analysis Dye Release Data River Reach.—Given the preceding analysis 
one can write for two stations, e.g., and 


and 


Numerical results are not shown for lack space, but results were comparable 
those Table except for the October 1977 data. 

For dispersion described diffusion equation, necessary for 
the motion tracer particle not correlated with its initial velocity 
(5). this basis initial convective period (5) length (18) defined. 
For the Mississippi River reach below St. Paul, the initial convective length 
about 1.9 miles considering only gravity-driven flow and wind effects. 
The measurements used herein were taken distance more than the convective 
length downstream from the injection point. Distances from the dye injection 


(6) (7) (9) (11) 
10.8 18.8 33.3 46.0 9.2 27.1 
0.00472 0.00198 0.000496 0.000331 0.00224 0.000817 
16,890 28,510 48,570 62,830 28,510 67,050 
(5,150) (8,690) (14,800) (19,150) (8,690) (20,430) 
447.1 325.1 362.9 505.8 1,020.4 
(34.77) (41.54) (30.20) (33.72) (47.00) (94.80) 
181.6 316.1 559.8 773.3 313.6 923.9 
(51.4) (89.5) (158.5) (219.0) (88.8) (261.6) 
11.3 7.14 6.26 4.33 6.00 1.84 
10.3 18.1 32.8 45.3 9.0 26.5 
0.081 0.081 0.081 0.081 0.081 0.081 
(0.265) (0.265) (0.265) (0.265) (0.265) (0.265) 
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CIS 
Model 


Concentration 


Time 
July 1976 


FIG. Versus Time Predicted CIS and Models; Dye Release 
Study, July 1976 


Model 
Model 


Concentration 


Time (days) 
May 


FIG. 8.—Concentration Versus Time Predicted and Models; Dye Release 
Study, May 1977 
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Model 
Model ----- 


Concentration 


Time (days) 
Oct. 


FIG. 9.—Concentration Versus Time Predicted CIS and Models; Dye Release 
Study, October 1977 


point station (Fig. and Tables and were 2.7 mile, 3.2 mile, and 


5.4 mile July 1976, May 1977, and October 1977, respectively mile 
1.609 km). 


Predicted C(t) Curves.—The parameters and used Eq. for 


the model predictions are those listed Table station station. The 
loss rate, was that given Table station station. Predictions based 
the model are very similar those made the CIS model shown 
Figs. and The predictive capabilities both models for the three 
case studies the Upper Mississippi are about equal. 


Riverine Versus Time Curves 


Skewness.—The CIS and models simulate concentrations versus time curves 
based bulk parameters derived estimated from observed concentrations, 
mean travel time and variance. Neither provides sufficient description 
the skewness the time concentration curves. 

The skewness the CIS and models directly related variance (see 
and and 28) and cannot adjusted independently. This deficiency 
more serious the observed time concentration curves become more skewed, 
seen the May 1977 data versus the July 1976 data for the Mississippi 
River study. The model results may very well agree with each other but differ 
from actual measurements when the distribution very skewed. The agreements 
and discrepancies shown Figs. and are directly related the magnitude 
the skewness coefficients shown Tables and 

strongly skewed distribution the time lag from the time maximum 
concentration the mean travel time (centroid) contributes heavily the 
magnitude the variance. Both the CIS and the model impose relatively 
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small skewness the predicted time concentration curves. When model variance 
matched with measured variance very skewed distribution, the results 
are shown Fig. (1) Lower predicted peaks the time concentration 
curves; and (2) wider spread the time scale. remedy this deficiency 
additional parameter (degree-of-freedom) would have introduced both 
models. 

The physical reasons for the skewness riverine time concentration curve 
are not the subject this paper, but few observations relative the Mississippi 
River study will follow. 

Concept Effective Volume Riverine Mass-Transport.—First, noted 
that the mean travel time, for the July 1976 and the May 1977 survey are 
not inversely proportional the river flow rate they would expected 
be, the travel time coincided with the hydraulic residence time (river volumes 
being essentially constant). the May 1977 survey, the dye cloud travels 
considerably faster than would justified the flow rate. Sampling location 
river cross section may also affect the concentration versus time curve 
but not the degree shown the data. 

The volumes, and kV, resulting from the application the CIS model 
the data and shown Table vary considerably from the July 1976 
the May 1977 survey. This direct disagreement with the real river volumes, 
which according stage-discharge relationships are nearly identical. One 
therefore led consider some sort short-circuiting effect the system, 
which can described introduction the The effective 
volume the volume obtained from the field dye experiment application 
the CIS model. measure that portion river system through 
which the dye travels. will normally less than the real volume the 
river system. 

the Mississippi River reach used this case study, vertical temperature 
stratification during the daytime has been observed under low flow (20). This 
daytime stratification accelerates transport the surface layer where the dye 
was introduced and measured. (All water samples for dye concentration 
measurements were retrieved within less than in. from the water surface.) 

believed that intermittent (daytime) stratification also responsible for 
the unusually high skewness the observed May 1977 time concentration curves. 
Stratification allows transport one layer while material left behind the 
other layers. 

Estimation Parameter from River dye tracings are not 
available cannot obtained for river reach, estimate the parameter, 
the CIS model can made present from river geometry and flow 
parameters follows. Equating the relationships, Eqs. and 27, and substituting 


The solution this quadratic equation for depends the numerical values 
the two dimensionless numbers, V,/Q and D/Ux. The former can 
calculated knowing the volume the river reach modeled, the kinetic 
rate coefficient the material, and the flow rate. The numerical value the 


HY6 CELLS-IN-SERIES SIMULATION 693 


Peclet numer, Ux/D, can found from river characteristics using any the 
empirical semi-empirical relationships for proposed other investigators 

For example, Liu and Cheng (16) proposed the relationship 
which cross-sectional area V,/x; mean depth the river reach; 
and coefficient determined from field data and ranges from 
according Liu and Cheng (16), with the order 0.6. Then, 
gravity; hydraulic radius the river, approximately equal the depth 
for wide rivers; slope the energy gradeline; and the width 
the river the reach under consideration, The preceding equation 
holds only dispersion dominated lateral variation velocity and not 
vertical velocity gradients. Substitution gives 
ghS 
which can used Eq. compute value for simplified solution 


Further studies the relationships between and river characteristics are 
needed. 


Riverine transport dissolved suspended material composite 
advection along river thalweg, transverse mixing towards the river banks, 
entrapment backwaters, and other separated flow regions created shoreline 
and bottom irregularities and other not well-identified processes. The prediction 
and analysis riverine transport for water quality modeling has been preoccupied 
with the application the equation, probably because early on, Taylor 
(23) showed how the equation can used analyze mass transport 
pipes. Elder (4) expanded its applicability wide and straight open channels 
and others later added transverse mixing (7,9,28) entrapment and dead zones 
(10,19,27), bends (6,8), etc. The geometry natural river often complicated. 
Spatial and time variations velocities and turbulence natural riverine flow 
systems are usually very pronounced. Simple one-dimensional models have 
been applied riverine systems incorporating the complicated aforementioned 
processes only one parameter, the longitudinal dispersion coefficient. Other 
models using expanded forms the equation lump parameters other 
than dispersion coefficient are worth consideration describe natural systems. 
this end, mixed CIS model has been explored this paper. 

The lump parameter the CIS model the number cells into which 
river reach subdivided. Each cell well-mixed. Dissolved material and 
water cascade from cells given river reach was kept constant the case 
study described herein but can made function flow rate where required. 
For water quality modeling, worth noting that biological chemical kinetics 
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are easily incorporated into CIS model. The mathematical form the governing 
equation, Eq. which ordinary linear differential equation with constant 
coefficients, makes possible express the response any function, 
means the convolution with the response unit impulse (21). 
Physically, this means that concentration cloud arbitrary shape can 
traced downstream superposition number unit impulses simulating 
the concentration cloud. Such application becomes very important the 
numerical analysis hydrologic transport branching systems, such shallow 
river impoundments, estuarine networks, ship canals, etc. 

has been shown this study that CIS model provides good description 
riverine transport one-dimensional model. Either model reproduces 
bulk features the mass transport (mean travel time and variance) without 
necessarily particular success the prediction instantaneous concentrations. 
The bulk parameters are the total mass dissolved suspended material, 
the mean travel time, the variance and the skewness, respectively, the time 
concentration curves. The definitions these four parameters have been given 
Eqs. and explicit forms applicable the CIS and the models, 
Eqs. 13-16 and 25-28, respectively. 

time-variable flow rates are used water quality model, the CIS 
model can easily expanded handle that feature. 

Numerical solutions the equation have similarity results from 
CIS model, because both cases river reach subdivided into finite, well-mixed 
volumes, and concentrations are obtained discrete values. Dispersion (numeri- 
cal dispersion) then becomes undesirable by-product the solution. CIS 
models, the dispersion due discretizing concentrations used simulate 
real dispersion. 

CIS models produce skewness without additional terms the transport 
equation. Retention material the cells integral part the model 
formulation. Skewness directly related the number cells used. 

The CIS model does not use diffusion-type term describe the dispersion 
process which introduces second-order derivative the conservation mass 
equation. The dispersion process related simple parameters such number 
cells, effective cell volume, and volumetric flow rate the CIS model. 

The actual volumes river reach are not required for the application 
the CID model. Instead, the necessary parameters, i.e., number cells and 
effective volume, can derived from dye-tracing experiment. fact, the 
use actual volumes based river geometry instead effective volumes 
may lead erroneous results travel times and longitudinal spread. 

dye tracings have not cannot performed river reach, estimates 
the parameter required for CIS model can made from river parameters 
using Eq. and any the empirical semi-empirical relationships proposed 
the literature for (e.g., Refs. 15, 16, 17, 19). Results such given 
Eq. can used until further studies have been conducted. 


The CIS models may most useful for slowly moving and highly dispersive 
rivers. 
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The following symbols are used this paper: 
area under concentration versus time curve (kilograms second 


per cubic meter); 
cross-sectional area (in square meters); 


width (in meters); 
concentration dissolved (tracer) suspended material (in 
kilograms second per cubic meter); 
Maximum value (in kilograms per cubic meter); 


dispersion coefficient (in square meters per sec- 
ond); 

dimensionless parameter, function Peclet number; 

mean depth (in meters); 

(in kilograms per cubic meters per second); 
rate material removal (in seconds); 

number cells river reach (—); 

total dissolved (or suspended) mass (in kilograms); 
released mass (in kilograms); 

Peclet number Ux/D; 
volumetric river flow rate (in cubic meters per second); 


time since material release (in hours); 

time centroid C(t) curve (in hours); 
mean river velocity (in meters per second); 
mixed cell volume (in cubic meters); 
material input into cell (in kilograms per second); 
distance from point material release (in meters); 


Subscripts 
1-2, 1-4, ... 
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D/x’; 

(Q/V) 

coefficient; 

Dirac delta function; 

skewness C(t) function (in cubic hours); and 
standard deviation C(t) function (in hours). 


refer river station number; and 
refer river reach between stations. 
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MINIMUM SPECIFIC ENERGY COMPOUND 
CHANNEL 


INTRODUCTION 


Analysis open channel flow the application the energy principle 
often clarified and aided the concept specific energy, which was introduced 
Bakhmeteff (1) 1912. Specific energy defined for one-dimensional 
open-channel flow the height the energy grade line above the channel 
bottom. leads classification open-channel flow into subcritical and 
supercritical flow regimes, distinguished flow depths that are respectively 
greater less than the depth which specific energy minimum (critical 
depth). mathematical consideration minimum specific energy gives rise 
the definition Froude number having value unity critical depth. 
The value the Froude number greater than unity for supercritical flow 
and less than unity for subcritical flow. 

The occurrence critical depth and its associated minimum specific energy 
considerable practical importance hydraulic engineers. one type 
channel control which may provide the boundary condition for computation 
water-surface profiles steady, gradually varied flow. Water-surface profile 
computations are integral part water resources investigations involving 
flood-plain delineations, evaluation flood control measures, and the design 
irrigation and drainage channels. 

Petryk and Grant (9) show that the determination critical depth channels 
with overbank flood-plain flow (compound channels) can troublesome. 
Customary definitions the Froude number generally not indicate critical 
depth the point minimum specific energy. addition, there are some 
compound-channel geometries which produce specific-energy diagrams with two 
points minimum specific energy. the purpose this paper present 
analytical formulation compound-channel Froude number which correctly 
identifies the occurrence points minimum specific energy for flow 
compound open channels. The proposed compound-channel Froude number can 
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used conjunction with existing computer programs for water-surface profile 
computation (5,13,16), and necessarily limited the same simplifying assump- 
tions that are associated with the conventionally used, one-dimensional equation 
steady, gradually varied flow (17). The results experimental investigation 
laboratory flume are also presented, demonstrating the existence two 


points minimum specific energy and identifying these points the proposed 
compound-channel Froude number. 


For simple channel nonrectangular section and uniform cross-sectional 
velocity distribution, the Froude number defined 


which discharge; the top width the water surface; acceleration 
gravity; and the cross-sectional area flow. For compound channel 
customary include the kinetic energy flux correction coefficient, 
the definition specific energy. result, appears follows the 
definition the Froude number assuming constant with depth: 


For natural channels with overbank flow, often assumed that the major 
contribution the large difference mean velocity between main channel 
and overbank sections. comparison the nonuniformity the velocity distribu- 
tion within each subsection can neglected. 

Two major problems arise the computation one-dimensional, steady, 
gradually varied flow profiles compound channels result using the 
Froude numbers First, incorrect solutions are generated when numerical 
methods are used solve the gradually varied flow equation written form 
involving the Froude number Second, incorrect solutions may accepted 
when the standard step method used compute water-surface profiles near 
critical depth. These difficulties are the result neglecting the variation 
with depth compound-channel flows. 

Consider the equation gradually varied flow the following form: 


which the rate change depth flow with respect distance 
along the channel; the bed slope the channel; and the slope 
the energy grade line. Prasad (10) has proposed numerical solution procedure 
for Eq. which can applied natural channels. addition the assumption 
that constant, the assumptions involved obtaining Eq. include: 
lateral flow, hydrostatic pressure distribution, constant bed slope, and 
straight, very wide channel, alternatively, approximately prismatic channel 
(17). Because the variation with depth and thus with distance along the 
channel has been neglected, application Eq. gradually varied flow 
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compound channel will lead incorrect water-surface elevations. The 
denominator the term the right-hand side Eq. arises from consideration 
the variation specific energy with depth, portion which due 
changes with depth compound-channel flow. Furthermore, the use 
can cause the right-hand side Eq. become indefinite depth 
that does not correspond the actual critical depth. 

alternative Eq. water-surface profiles are computed natural 
channels the standard step method (6) which the specific energy computed 
explicitly. this case, does not appear the equation solved, but 
used instead indicate whether the solution the supercritical subcritical 
flow regime. For compound channels, neither nor correctly indicates the 
flow regime. Thus, incorrect solutions the energy equation can accepted 
when the depth near critical depth. 


Previous Investigations.—Previous investigations the problems associated 
with defining the Froude number compound channels are limited. Numerous 
laboratory investigations compound-channel flow have been undertaken 
(8,11,15), but the focus these experiments has been the quantification 
changes the boundary shear stress distribution resulting from momentum 
exchange between the main channel and floodplain. The Federal agencies which 
maintain and use water-surface profile programs recognize the Froude-number 
difficulties compound channels described the previous section this 
paper, and they examine these difficulties their user’s manuals. The Soil 
Conservation Service (16), e.g., warns differences much between 
the critical depth determined (Eq. and the critical depth determined 
minimum specific energy. 

The Corp Engineers (5) presents algorithm solve for the depth 
corresponding minimum specific energy when their water-surface profile 
program attempts obtain solution close critical depth. The depth 
minimum specific energy compared with the profile depth check the flow 
regime rather than using the Froude number check. 

The United States Geological Survey (USGS) (12) proposes the use 
index Froude number based the Froude number the subsection carrying 
the greatest discharge. The index Froude number thought the USGS 
better reflect the flow regime the entire cross section, but also recognized 
having limitations. The USGS does not consider the index Froude number 
true Froude number, but rather warning flag that identifies possible 
flow-regime problems. later version the USGS Water Surface Profile Program 
incorporates routine determine the depth minimum specific energy. 

Petryk and Grant (9) have proposed discharge-weighted Froude number 
without experimental corroboration order eliminate the computational 
problems associated with the occurrence two points minimum specific 
energy compound-channel flows. Although their proposed Froude number 
succeeds doing this identifying only one value critical depth, 
nevertheless somewhat arbitrary and divorced from the concept minimum 
specific energy. 

Clearly, the Froude number should formulated reflect the specific energy 
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curve under consideration and should indicate critical depth the point (or 
points) minimum specific energy. Such Froude number would produce 
correct numerical solutions the gradually varied flow equation (Eq. and 
would eliminate the need for time-consuming routines used solve for the 
depth minimum specific energy standard step water-surface profile 
computations. 


Derivation and Formulation.—The specific energy, for one-dimensional 
compound-channel flow given 


which the depth flow. The kinetic energy flux correction coefficient, 
defined 


which the velocity through the element area, dA; and the mean 
cross-sectional velocity (3,6). Alpha thus measure the nonuniformity 
the velocity distribution. For computational purposes, flow conventionally 
divided into channel and overbank subsections appropriately located vertical 
lines which are assumed not transmit shear stress from one section flow 
another, and which not contribute wetted perimeter. Wright and Carstens 
(15) have suggested that the wetted perimeter the subsection dividing line 
retained for the main channel, and that the shear stress applied the 
main-channel flow section the overbank section considered. Regardless 
the manner which the main flow-flood-plain interaction treated, the 
basic assumption the computation previously mentioned, that 
the contribution the nonuniformity the velocity distribution within each 
subsection negligible comparison the variation mean velocity between 
subsections. Eq. applied with this assumption compound channel 


which has been divided into subsections, the kinetic energy flux correction 
coefficient becomes 


which the conveyance the ith subsection; the area the ith 
subsection; and the conveyance the total cross section (3,6). 
The subsection conveyance computed from the Manning equation follows: 


which a,/p,) the subsection hydraulic radius; the subsection 
wetted perimeter; and the subsection value. the system units, 
the constant replaced unity. 


1.49 
n, 


COMPOUND OPEN CHANNEL 703 


The point (or points) minimum specific energy obtained differentiating 
Eq. with respect and setting the derivative equal zero. Because both 
and area are functions depth, the differentiation produces (14) 


Noting that /dy and that rearranging terms, the following expression 
obtained: 


The left-hand side Eq. unity the point minimum specific energy; 
therefore, compound-channel Froude number can defined from Eq. 


1/2 


the point minimum specific energy will have value 


FIG. 1.—Definition Sketch for Evaluation 


With the exception da/dy, all the terms the right-hand side 
Eq. are routinely determined water-surface profile computations. Evaluation 
da/dy can achieved differentiating Eq. with respect shown 
Appendix the derivative becomes 
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Eqs. 12-14, the top width the ith subsection; and dp,/dy the 
rate change wetted perimeter with respect depth flow the ith 
subsection. Evaluation simplified the fact that the cross-section 
geometry natural channels defined ground points connected with straight 
lines. The definition sketch Fig. (which portion right overbank 
subsection) shows the water surface intersecting the line segment de. This line 
segment makes contribution the subsection wetted perimeter. The 


rate change wetted perimeter with respect depth constant along 
de, and therefore can evaluated 


The terms and are generally determined when computing the geometric 
properties cross section for use water-surface profile program. 
should noted that the water surface point dp, should evaluated 
for the line segment de, but the water surface point should 
evaluated for the line segment cd. situations where the water surface 


FIG. Cross Sections for Evaluation Specific Energy and Froude 
Numbers: (a) Cross Section (b) Cross Section 0.3 


does not intersect the wetted perimeter subsection (e.g., the boundary 
between the main channel and overbank section above bankfull stage), dp, 
zero. For subsection where the water surface intersects both left and 
right bank (e.g., the main channel below bankfull stage), dp,/dy the sum 
Ap/Ay for each the banks. 

The working equation for the compound-channel Froude number can 
obtained substituting Eq. into Eq. and simplifying: 


the Manning’s value considered vary with depth flow any 
subsection, and can written reflect the variation: 


dp, 


6ft 
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Evaluation.—The behavior the compound-channel Froude number, may 
evaluated examining the specific-energy diagrams two idealized, 
symmetric cross sections, each conveying 5,000 cfs (142 Cross section 
[Fig. from Petryk and Grant (9). Fig. the specific-energy curve 
for this cross section reveals two points minimum specific energy depths 
flow approx 6.8 (2.07 and 5.3 (1.62 m). These points are indicated 
and C2, respectively, Fig. 

(Eq. 16) for this cross section plotted Fig. along with (Eq. 
and (Eq. 2). expected, all three equations produce the same curve below 
top bank (simple channel situation), but only Eq. for correctly locates 
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FIG. 3.—Specific Energy for Cross Section Conveying 5,000 cfs cfs 0.028 


Cl, the upper depth minimum specific energy (6.8 2.07 m), and connects 
with the lower curve the top bank depth. 

The shape the Froude number curve independent the discharge, and 
the fiducial point can shifted left right varying the discharge. 
This means that once plotted for particular cross section and discharge, 
points minimum specific energy for other discharges may determined 
without the necessity constructing new specific-energy diagrams. effect, 
the variable provides universal horizontal scale for Fig. which depends 
only the conveyance and geometric properties the particular cross section. 
Thus, for given depth flow, the critical discharge, can computed 
taking the reciprocal the corresponding value because 
for the given depth equals for the critical condition. 

Cross section presented Fig. 2(b) and differs from cross section 
only that the flood plains have 100:1 slope toward the channel. The 


HY6 705 


706 JUNE 1981 HY6 


specific-energy diagram cross section (Fig. reveals single point 
minimum specific energy below top bank the same depth flow for 
cross section (point C2). The three Froude number curves shown Fig. 
for cross section are again identical below top bank, but (Eq. 
and (Eq. each indicate another point minimum specific energy above 
top bank depths flow 6.5 (1.98 and 6.8 (2.07 m), respectively. 
The occurrence these false points minimum specific energy more 
serious deficiency Eqs. and than the errors critical depth shown 
Fig. 

evident from these two examples that the Froude numbers generated 
and are not acceptable for use the gradually varied flow equation 
(Eq. 3). Neither definition Froude number faithfully reflects the specific-energy 
diagram overbank flow situations, and either would produce divergence from 
profile solution. equally evident that Eqs. and are not satisfactory 
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FIG. 4.—Froude Numbers for Cross Section Conveying 5,000 cfs cfs 0.028 


for checking the flow regime the standard step method. Only (Eq. 17) 
accurately reflects the specific-energy diagram and indicates the correct flow 
regime. The experimental investigation into the occurrence two points 
minimum specific energy the following portion this paper offers guidance 
for the interpretation the flow regime between the two points minimum 
specific energy, and C2, cross section (Fig. 3). 


EXPERIMENTAL INVESTIGATION 


The experimental investigation consisted measuring point velocities 
compound-channel cross section which was formed constructing single 
rectangular overbank section laboratory flume. Sufficient point velocity 
measurements were made eight different depths flow (at approximately 
the same discharge for each depth) compute the discharge, mean velocity, 
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kinetic energy flux correction coefficient, and specific energy for each depth. 
Complete details the experimental procedure are given the first writer 
(2). 

The experiments were conducted tilting steel flume (24.38 long, 
3.5 (1.07 wide, and 1.5 (0.46 deep. This flume was also used 
Tracy and Lester (14) and details its construction are given them. The 
overbank section was constructed 3/4-in. exterior plywood and two-by-six 
fir framing lumber, resulting the channel dimensions shown Fig. All 
wood was coated with sanding sealer and exterior acrylic-latex paint. The overbank 
section was attached the flume with silicon adhesive. 

Point velocities were measured with 0.072-in. (1.83-mm) outside diameter 
pitot-static tube operated conjunction with differential pressure transducer. 
Data collection, reduction, and analysis were accomplished with HP9825A 
desktop computer controlling digital voltmeter which measured the voltage 
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FIG. 5.—Specific Energy for Cross Section Conveying 5,000 cfs cfs 0.028 


output from the pressure transducer and preamplifier. Point velocity measure- 
ments were made station (19.81 downstream the flume entrance. 
Preliminary measurements were made station (18.29 downstream. 
Comparison dimensionless profiles velocity between the two stations 
indicated that the flow was fully developed. 

The preliminary experiments indicated that discharge 1.7 cfs (0.048 
would produce specific-energy curve with two points minimum specific 
energy. estimate the error setting the discharge 1.7 cfs (0.048 
included the calibration error the Venturi meter used measure the discharge 
and also included estimate the error introduced observed fluctuations 
the Venturi-meter manometer during the course experimental run. The 
estimated error discharge was the order +3%, which was the same 
range error observed between individual discharges determined from the 
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Ver‘uri meter and the discharges determined integration the point velocity 
measurements. 

Establishing truly uniform flow profile for the experimental runs proved 
impossible. Any discharge flowing near the depth corresponding minimum 
specific energy, these were, could expected inherently unstable. 
The instability was exacerbated the variations the overbank surface, which 
were the order +0.01 (0.3 cm). Standing waves the overbank section 
and cross-hatched water surface the channel thwarted efforts achieve 
uniform water-surface profile. result, the adopted experimental procedure 
was establish profile close uniform possible such that the desired 
depth flow was obtained where the point velocities were measured. 
The maximum observed change depth for overbank-flow runs was approx 
0.05 (1.5 cm) between the channel entrance and the measuring station where 
the flow depth was 0.567 (17.3 cm). For larger depths flow, the water-surface 


- 
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° 

— 
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Froude Number 


FIG. 6.—Froude Numbers for Cross Section Conveying 5,000 cfs cfs 0.028 


profiles tended more stable and more nearly uniform. profile 
depth flow 0.7 (21.3 cm) was established demonstrate that uniform 
profile could obtained the downstream reach the flume the depth 
flow was sufficiently greater than the depth corresponding minimum specific 
energy. 


Table presents the values area, discharge, kinetic energy flux correction 
coefficient, and specific energy computed from experimental measurements for 
each the eight reported runs. Runs and are not reported the table 
because operational difficulties during each run. apparent from the 
results presented Table that the depth increased for those experimental 
runs with overbank flow, the proportion the total discharge the overbank 
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section increased. should also noted that the values for the main 
channel alone are measurably larger than 1.0 because the narrowness 
the main channel section. 

Observations the water surface for the four experimental runs with overbank 
flow indicated greater instability the depth flow decreased. The water-surface 
instability was manifested standing waves the overbank section and 
choppy, cross-hatched water surface the channel section. Beginning the 
upper depth minimum specific energy (run and continuing with decreasing 
depth, the standing wave fronts the overbank section were perpendicular 
the mean flow direction and then were bent downstream into cross-hatched 
pattern the channel section characteristic supercritical flow. The surface 
instability continued increase for the experimental runs depth decreased 
below top bank. The fact that the water surface was unstable for experimental 
runs and the first two runs below top bank Table suggests that 
the upper point minimum specific energy could considered the limit 
subcritical flow for situations which two points minimum specific energy 
occur water-surface profile computations. 


Overbank Section 


3.490 


FIG. 7.—Cross Section Flume and Overbank Section, Looking Downstream 
0.3048 


The experimental specific-energy data Table are plotted Fig. 8(a). 
Although the variation discharge from run run causes some scatter 
the plot, there evidence two points minimum specific energy. The 
experimental values plotted Fig. 8(b) show little scatter and indicate 
that primarily function depth flow. This observation suggests that 
specific-energy diagram for single value discharge can constructed 
substituting the average discharge eight runs (1.692 cfs 0.048 
into Eq. while using the experimental data for all other variables. Fig. 
presents the resulting average specific-energy diagram. The two points minimum 
specific energy are more clearly apparent this figure. 

The concept computing Froude number for the flow subsection 
compound channel has already been mentioned with regard the USGS 
index Froude number (12). The subsection Froude numbers (computed with 
Eqs. and for the experimental data this investigation are presented 
the index Froude number these experimental runs because the channel 
the subsection with the largest discharge. All four depths flow above top 
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Channel 
cubic 
square feet per 
feet second 
(5) (6) 
0.001018 
0.001128 
0.001485 
0.002096 
0.001903 
0.002118 
0.003300 
0.004455 


bank are subcritical based the index Froude number, but shown 
Fig. the two lower overbank depths are not subcritical. For this experimental 
investigation, the index Froude number does not correctly indicate the flow 
regime compound-channel flow. 

Petryk and Grant (9) apply the concept subsection Froude number 
obtain their weighted Froude number which given 


which the subsection discharge; and the subsection Froude number 
computed Eq. Values for the experimental data are presented 
Col. Table the case the index Froude number, the weighted 
Froude number does not correctly indicate the flow regime. 


The proposed compound-channel Froude number cannot directly determined 
from the experimental data. Attempts use Eq. fail because difficult 
determine da/dy from the limited number experimental data points. Eq. 
fails because the slope the energy grade line not precisely known, 
which means that the subsection resistance coefficient and thus the conveyance, 
k,, cannot determined from the experimental data. had been possible 
establish uniform flow condition for each run, the energy gradient would 
parallel the flume slope, and the conveyance for each subsection could 
computed from the experimental data alone. The compound-channel Froude 
number can only determined indirectly through independent prediction 
the experimental results. 


Working the same flume used the present investigation, Tracy and 


Run 
(1) (7) 
0.650 1.084 
0.625 1.083 
0.600 1.082 
0.567 1.088 
0.533 1.093 
0.500 1.087 
0.467 
0.433 1.100 
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Overbank 
cubic A,in cubic 
square feet per square feet per 
feet second feet second 


(8) (9) (11) (12) 


Lester (14) experimentally determined friction-factor relationship for smooth 
rectangular channels the form 


Vf) 


which the Darcy-Weisbach friction factor; and the Reynolds number. 
assumed that Eq. valid when applied independently each channel 


subsection, the friction factor, f,, can determined for the ith subsection. 
The mean velocity the ith subsection, v,, then given 


which the hydraulic radius the ith subsection; and the slope 
the energy grade line. Because the values and obtained from 
and must such that the subsection discharges sum the average 
measured discharge, Q,,, 1.692 cfs (0.048 the following equation must 
satisfied: 


= To ft 


has been implicitly assumed that the same for all subsections. 
20, 21, and can solved iteratively for the friction factor and velocity 
each subsection for given total discharge and depth. The iterative solution 
procedure given detail the first writer (2). 

The velocities, v,, were calculated the procedure just described for the 
mean measured discharge 1.692 cfs (0.048 was assumed that the 
imaginary vertical boundary between the main channel and overbank section 


a 
(13) 
0.294 0.411 1.108 0.927 1.774 1.192 
0.231 0.326 1.132 0.840 1.714 1.198 
0.169 0.230 1.169 0.753 1.726 1.224 
0.083 0.087 1.340 0.635 1.680 1.238 
0.519 1.648 1.093 
0.487 1.676 1.087 
0.455 1.645 1.096 
0.422 1.671 1.100 
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made contribution wetted perimeter. Furthermore, the friction factors 
determined for each subsection were converted Manning’s values because 
the formulation for the compound Froude number, terms The 
values obtained exhibited slight variation with depth; however, facilitate 
the computations, constant values 0.009 and 0.010 were adopted for the 
channel and overbank sections, respectively. From the velocities and values 
for each subsection, the specific energy and compound Froude number were 
computed for series depths within the range measured depths. the 


TABLE 2.—Froude Numbers for Experimental Data 


(Eq. 
(4) 


feet 


e- channel 


Depth Flow, feet 


Depth Flow 


o-total 


068 072 


FIG. 8.—Specific Energy and Kinetic Energy Flux Correction Factor from Experimental 
Data 0.3 m): (a) Specific Energy, feet; (b) Alpha 


computation the specific energy and was assumed that each 
subsection had the value 1.0 rather than the measured value. this way, the 
computational procedure remained independent the measured data and was 
executed the same manner would expected when determining for 
natural river channel the course water-surface profile computation. 
The predicted specific-energy diagram shown Fig. 10(a), and two depths 
minimum specific energy are apparent, although each depth approximately 
2/100 smaller than the corresponding depths Fig. Fig. The 


Run feet (Eq. (Eq. 19) 
(1) (2) (3) (7) 
0.650 0.471 0.490 0.721 0.759 0.529 
0.625 0.508 0.529 0.821 0.873 0.586 
0.600 0.583 0.606 0.925 1.001 0.629 
0.567 0.675 0.704 1.017 1.177 0.692 
055 
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entire specific-energy curve Fig. 10(a) skewed slightly downward and 
the left when compared with the measured curve Fig. the average 
curve Fig. The predicted compound-channel Froude number curve Fig. 
exhibits the behavior typical for two points minimum specific energy, 


a 


top bank 


Depth Flow, feet 


+ 


O72 
Specific Energy, feet 


FIG. 9.—Experimental Specific Energy Curve for Average Discharge 1.692 cfs 


ao 


top bank top bank 


9 


Depth Flow, feet 


Depth Flow 


FIG. Predicted Specific Energy Experimental Flume for 1.692 cfs; (b) 
Compound Channel Froude Number for Fig. 10(a) cfs 0.028 0.3 


and correspondence with the predicted specific-energy curve expected. 
investigate the role that neglecting the transfer linear momentum 

the overbank section plays the skew the predicted specific-energy curve, 

the correction suggested Wright and Carstens (15) was considered. Although 
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055 
050 
045 045 
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the correction improved the agreement between the measured and computed 
discharges the overbank section, especially the larger depths, the effect 
the computed specific-energy curve was minimal because the relatively 
small changes which resulted from the correction. 

The skew the specific-energy curve most pronounced below top bank 
depth where transfer linear momentum the overbank does not occur. 
The skew this portion the curve can attributed selecting subsection 
values unity computing specific energy. should noted that the 
depths flow the flume were small compared depths flow normally 
found field situations. For this reason, the velocity head the flume makes 
large relative contribution specific energy, and any adjustment velocity 
head (such subsection has far more effect specific energy the flume 
than would the field. 

The same analysis can applied subcritical and supercritical flow regimes 
field situations where kinetic energy correction coefficients can much 
1.4 more the main channel (7). For subcritical flow where the velocity 
head small, a-adjustment velocity head would insignificant. For 
supercritical flow, the velocity head can 50% more the depth, and 
a-adjustment velocity head would have significant effect specific 
energy. This reasoning explains the increasing leftward shift Fig. 10(a) 
the depth flow decreases, and the implication that predicted specific energies 
and Froude numbers field channels under subcritical flow conditions would 
closer measured values. 


Existing formulations the Froude number (Eqs. and not accurately 
reflect the specific-energy curve for flow compound open channel and 
not correctly locate points minimum specific energy. compound-channel 
Froude number (Eq. 16) has been derived and has been shown accurately 
reflect the specific-energy curve flow compound open channel correctly 
locating points minimum specific energy. When applied simple channel 
with uniform velocity distribution, the compound channel Froude number 
identical Eq. the conventional definition Froude number. 

The compound-channel Froude number appropriate for use with the gradually 
varied flow equation (Eq. and provides the proper check the flow regime 
when used conjunction with the standard step method water-surface profile 
computation. The proposed Froude number subject the same assumptions 
that apply the equation gradually varied flow commonly employed 
water-surface profile computations. 

For some compound-channel geometries characterized wide, level flood 
plains, two points minimum specific energy can computed for certain 
discharges. Laboratory investigation one-dimensional flow demonstrates 
that this phenomenon can fact occur, and indicates that the upper point 


minimum specific energy may considered the proper limit subcritical 
flow. 
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Writing Eq. 


and diff with respect produces 


Noting that da,/dy t,, dA/dy and dK/dy the following 
obtained: 


a,/ 


Evaluate dk, /dy writing Eq. 


and differentiate with respect obtain 


Again noting that t,, and multiplying and dividing a,, the following 
obtained: 


Substituting Eq. into Eq. and simplifying, results Eq. 11. 
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The following symbols are used this paper: 


total cross-section area; 

subsection area; 

specific energy; 

Froude number; 

compound-channel Froude number; 

subsection Froude number; 

weighted Froude number; 

Froude number with kinetic energy flux correction; 
Darcy-Weisbach friction factor; 

subsection friction factors; 
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acceleration gravity; 

total cross-section conveyance; 
subsection conveyance; 
Manning’s value; 

subsection value; 

subsection wetted perimeter; 
total cross-section discharge; 
average measured discharge; 
subsection discharge; 

Reynolds number; 

subsection hydraulic radius; 
slope energy grade line; 

bed slope channel flume; 
total cross-section top width; 
subsection top width; 

total cross-section mean velocity; 
mean velocity associated with incremental area, dA; 
subsection mean velocity; 
distance along channel; 

depth flow; 

kinetic energy flux correction coefficient; 
increment wetted perimeter; 
increment depth; and 


subsection parameters compound-channel Froude number. 
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DIFFUSION- WAVE ROUTING 
CHANNEL NETWORKS 
Ali Osman and Ben Chie Yen,’ ASCE 


For subcritical flow open-channel network, mutual backwater effects 
exist among the channel branches joining junction. Therefore, the branches 
cannot treated individually when dynamic wave model (4,10) adopted 
route floods open-channel network. Ideally, the entire network should 
considered single unit and the flow all the channels and junctions 
should solved simultaneously. However, for network consisting more 
than few channels, this approach requires excessive amount computer 
storage and expense. 

the sequential type models (3,6) traditionally employed for routing floods 
through channel network, the difficulty caused the mutual backwater effect 
channels does not arise because the backwater effect from downstream 
ignored. The routing starts from the most upstream channels and carried 
towards downstream, channel channel sequence, satisfying the flow 
continuity requirement junctions. The downstream boundary condition 
channel required the computations treated somewhat casually since 
actually unknown the time seeking solution. However, Seviik and 
the second writer (8), and the writers (11) have shown that sequential models 
give unrealistic solutions when the downstream backwater effect significant. 

order reduce the high computer cost and storage requirements solving 
large network problems, successive decomposition technique called the over- 
lapping-segment scheme has been employed routing floods through channel 
networks (7,11). example given Fig. Solution sought separately 
and successively for each individual segment which consists one junction 
and the channels joining it. Consequently, the computer storage requirement 
substantially decreased. When the Saint Venant equations, together with the 
junction continuity and dynamic equations are employed describe the flow, 
the coefficient matrix involved the solution process not banded and therefore 
the matrix solution technique not most efficient. 

Prof. Civ. Engrg., Middle East Technical Univ., Ankara, Turkey. 

Civ. Engrg., Univ. Illinois Urbana-Champaign, Urbana, Ill. 61801. 

Note.—Discussion open until November 1981. extend the closing date one month, 
written request must filed with the Manager Technical and Professional Publications, 
ASCE. Manuscript was submitted for review for possible publication July 22, 1980. 
This paper part the Journal the Hydraulics Division, Proceedings the American 
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Another avenue reduce the computational costs flow routing use 
approximation the Saint-Venant equations (dynamic wave) instead 
the complete equations. well-known that only approximations the level 
diffusion wave model higher can account for the downstream backwater 
effect whereas the kinematic wave models cannot (10). Ponce, al. (5) 
investigated the conditions under which the diffusion wave can good 
approximation the dynamic wave model. 

this paper, nonlinear diffusion wave model simulate the unsteady 
flow dentritic network accounting for the downstream backwater effect 
presented. The overlapping-segment technique applied the diffusion wave 
equations which are written finite differences. The resulting set nonlinear 
algebraic equations for each segment solved using the Newton iteration 
method. The coefficient matrices obtained the solution process are banded 
and, therefore, the matrix equations are solved very efficient numerical 


FIG. 1.—Overlapping Segments 


technique. The results the proposed model compare favorably those 
complete dynamic wave and kinematic wave model. 


Gradually varied unsteady flow open channels mathematically described 
set one-dimensional shallow water equations commonly known the 


Saint-Venant equations. gravity-oriented coordinate system, these equations 
are written (9,10) 


which distance measured horizontally along the channel; flow 
water surface elevation measured vertically from horizontal datum; 
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momentum correction coefficient; gravitational acceleration; friction 
slope; and time. For practical purposes, usually taken unity. The 
friction slope for turbulent flow can estimated using Manning’s formula 


which the Manning roughness factor; hydraulic radius; and the 
constant for units and 1.486 for U.S. Customary units. Alternatively, 
other friction formulas such Darcy-Weisbach’s Chezy’s may also used. 

Hydraulic conditions junction may described the continuity equation 


and the dynamic equation 


2 v2 


which the storage within the junction; head loss through the junction; 
and flow velocity. The subscript stands for any one the in-flowing 
channels and represents the out-flowing channel. Other symbols are 
previously defined. 


Junctions have small storage volumes most open-channel networks for 
which the term ds/dt Eq. negligible. Thus, one can write 


Also, when the flows all the branches joining junction are subcritical, 
Eq. can approximated kinematic compatibility condition 


Diffusion Wave the diffusion wave approximation the Saint- 
Venant equations, the local and convective acceleration terms the momentum 
equation (i.e., the first two terms Eq. are neglected (10). Thus, Eq. 
simplified 


which may account for flows both positive and negative directions. 


Combining and yields 
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Finite Difference Equations.—Consider Y-segment shown Fig. 
example. the Y-junction, and are the two incoming channels and 
the outflowing channel. These channels are subdivided and 
cross sections, respectively, for the branches and into 
and computational reaches. Numerical solutions are sought 
the flow cross sections. Although the reach length, Ax, allowed vary 
the model, treated here constant for channel for simplicity. 


Spoce,x 
FIG. 3.—Computational Grid 


Inflow discharge hydrographs the upstream end channels and are 
given the upstream boundary conditions. other words, and are 
prescribed functions time where the first subscript denotes the channel and 
the second subscript stands for the flow section considered. downstream 
boundary condition for channel must also known. This condition can 
prescribed stage hydrograph stage-discharge relationship. 

With reference Fig. Eq. written finite difference quotients for 


Junction 
FIG. 2.—Discretization Y-Segment 
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j j j+l j+i 
A’giss + — — Bi+t Bui 


2At 


Stage which solutions are sought; and the previous time stage which 
solutions are already known. The time increment between the two stages 
denoted Ar. The unknown terms Eq. are those containing the superscript 
except for which specified the upstream boundary condition. 
The discharge, interior flow section computed using Eq. which 
written finite difference form 


j+i j+i 
iss h’s 
Ait! Ri?! 2/3 Ax 
( j+i j+i 


Since the flow area and the hydraulic radius are known functions the 
water surface elevation substituting Eq. into Eq. for ..., 


Similarly, for channel set equations with unknowns are 
obtained 


For channel set finite difference equations with unknowns 
for the interior flow sections 


The downstream boundary condition for channel expressed 


unknowns. Additional three equations are obtained from the junction 
conditions. For the Y-junction being considered, Eq. becomes 


Also, Eq. written separately for branches and yields 


h’s 


Making use Eq. and the cross-sectional geometry the channels, Eqs. 
16-18 are expressed 
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Solution Technique for Finite Difference Equations.—The Newton iteration 
method employed for the solution the simultaneous nonlinear finite difference 
equations. The procedure can summarized follows: 


set trial values, are assigned the unknowns 

Substituting into Eqs. 12-15 and 19-21, values are computed. 

Partial derivatives the functions appearing Eqs. 12-15 and Egs. 
19-21 with respect h’*' are taken and evaluated usingh’*' 

The matrix equation shown Fig. obtained using the values 
and af/ah’*'. The superscripts are not shown Fig. for the sake 
clarity. should noted that the variables, and Eq. 19, 
and the variable, Eqs. and are treated constants assuming 
their values obtained the previous iteration cycle while computing the 
partial derivatives the Only with this approxima- 
tion can banded coefficient matrix obtained Fig. 

The matrix equation Fig. solved for the corrections Ah. the 
magnitudes these corrections are within tolerable limits for all the flow sections, 


FIG. 5.—Hypothetical Network Used Examples and 


the trial values are accepted. Otherwise, the trial values are improved 
for all the flow sections, and new iteration cycle initiated. 


should noted that the coefficient matrix Fig. banded with 
most three nonzero elements the diagonal. This allows very efficient solution 
procedure which has been described elsewhere (1) detail and will not 
repeated here. 

Overlapping-Segment flow channel network containing 
more than one junction simulated, the proposed model employs 
overlapping-segment technique mentioned the shown 
each containing only one junction. The outflowing channel segment 
also in-flowing channel the immediate following segment. For instance, 
Fig. the out-flowing channel segment serves also the inflowing 
channel IIIC for segment III. 

Each segment treated separate unit the computation. Computations 
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start from upstream segments and progress towards downstream. For example, 
Fig. first, segments and are considered separately with the downstream 
boundary conditions for channels and IID assumed approaching the 
normal flow for long channel. Finite difference equations for each these 
two segments are obtained and solved respectively manner similar those 
previously described for the simple Y-segment Fig. The only exception 
that Eq. replaced uniform flow equation corresponding the 
assumed downstream boundary condition. Solutions obtained are assumed valid 
for channels IB, IC, IIB, and IIC. However, solutions for flow channels 
and IID are discarded and then recomputed channels and IIIB 
segment III. Naturally, the algebraic sum outflows from channels and 
each time level constitutes the upstream inflow for channel (IIIC). 
Similarly, the sum outflows from IIB and IIC gives the upstream boundary 
condition for IID (IIIB). Knowing the upstream inflows for channels IIIC and 
IIIB, computations are performed for segment III. 


TABLE 1.—Properties Example Open-Channel Network 


Length, Width, 
Channel feet feet Manning 
number (meters) (meters) 


n 
(1) (2) (4) (5) 


The single overlapping scheme just described may not sufficiently accurate 
the downstream backwater effect serious and significant backwater propagates 


beyond one upstream channel. this case double overlapping scheme may 
adopted. 


Comparison with 


order test the validity the proposed diffusion wave model for simulation 
flow open-channel networks, comparisons with two other methods, namely, 
dynamic wave and kinematic wave models, have been made. The two models 
are briefly described follows. 

Implicit Nonlinear Dynamic Wave four-point fully implicit scheme 
which adopted this study for dynamic-wave routing modified version 
that proposed Baltzer and Lai (2) incorporate the simultaneous solutions 


1968.5 0.0005 16.4 0.0138 
(600) (5) 
1968.5 0.0005 16.4 0.0207 
(600) (5) 
1968.5 0.0005 16.4 0.0207 
(600) (5) 
1968.5 0.0005 16.4 0.0138 
(600) (5) 
1968.5 0.0010 26.2 0.0141 
(600) (8) 
1968.5 0.0010 32.8 0.0125 
(600) (10) 
7 
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the Saint-Venant and the junction equations. this scheme, the finite difference 
quotients for reach between flow sections and are written 


l j+ + 
Fisos > (Fi ‘+ 


i+0.5 t 


which any reach function. The set finite difference equations obtained 


100 
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FIG. 6.—Upstream Inflow Hydrographs for Example 
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FIG. 7.—Upstream Inflow Hydrographs for Example 


substituting Eqs. 22-24 into Eqs. and contains two unknowns, namely 
and for each flow section along channel. When simple channel 
network composed three junctions and seven branches shown Fig. 
considered, each time step, the finite difference forms and 
provide 14(J nonlinear equations with unknowns the number 


400 
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200 
0 Channets 
300 
Cha 
100 
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flow sections selected along each channel The additional equations 
necessary for solving the problem are provided four upstream boundary 
conditions for channels IB, IC, IIB, and one downstream boundary condition 
for channel IIID, and nine junction equations, namely, Eqs. and applied 
junctions II, and III. The solution obtained using the Newton iteration 
technique. The coefficient matrix obtained each iteration cycle has rank 
14J, and not banded. This makes the solution difficult and costly 
for small networks like that Fig. and unmanageable for larger systems. 

Nonlinear Kinematic Wave Model.—In the kinematic-wave approximation, the 
inertia well pressure terms the momentum equation (Eq. are neglected 
and thus, the friction slope, S,, approximated the bottom slope, S,, 
the channel (10). The downstream backwater effects channels are neglected 
and each channel treated individually. Computations start from the most 
upstream channels for which upstream boundary conditions are specified, and 
progress towards downstream channel channel cascading manner, satisfying 
the flow continuity requirement junctions. Along channel reach with 
Eqs. and are combined yield 


For channel along which flow sections are considered, applying Eqs. 
22, 23, and Eq. and writing for grids, together with upstream 
boundary condition, set nonlinear algebraic equations with unknowns 


obtained. The unknowns are for ..., This set nonlinear 
equations are solved using the Newton iteration method. The coefficient 
matrices involved are banded and thus, the solutions are obtained very 
efficient manner. 

Comparison Models.—Two hypothetical examples are presented here 
demonstrate the validity and efficiency the proposed diffusion wave model 
compared the other two models. The hypothetical network used the 
examples consists two junctions and six branches shown Fig. All 
branches are assumed rectangular cross section. The dimensions 
the channels are listed Table The hypothetical floods that are routed through 
the network the two examples are shown Figs. and respectively. 
The initial condition for both examples steady flow condition corresponding 
channel Downstream boundary condition the exit channel specified 
uniform flow equation assuming this channel hydraulically long. time 
increment sec selected both examples for all the three models. 
constant reach length (196.9 ft) adopted for the diffusion 
and dynamic-wave models. The reach length reduced (98.4 
ft) the kinematic-wave model since numerical convergence could not 
achieved with (196.9 ft). 

The computed discharge hydrographs using the three models are shown 
Figs. and respectively, for Examples and II. can observed from 
these figures, the results the proposed diffusion-wave model are good 
agreement with those the dynamic wave model. Conversely, the kinematic 
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wave solutions are not reliable. The outflow hydrographs computed for channels 
and Example and shown Fig. clearly demonstrate the downstream 
backwater effects. Example constant upstream inflow (71 
cfs) equal the base flow rate adopted for channels and shown 
Fig. Thus, steady flow condition would prevail both channels and 
the time variant backwater effects the other channels the network 
did not exist. Indeed, the kinematic wave method which ignores the backwater 
effects from downstream predicts steady flow these two branches shown 


Dischorge , 


FIG. 8.—Comparison Results Example (a) Channel Channels and 
(c) Channel and (d) Channel 


Fig. However, actuality, the flood wave travelling through channel 
raises the water surface the junction where branches and join. 
This decreases the hydraulic gradient channels and causing discharges 
lower than the constant upstream inflow and base flow rate. Naturally, from 
continuity requirements, the channel storage increased during the period 
low discharges. the downstream backwater recedes with time, the excess 
water released from the channel storage, and discharges higher than the constant 
inflow rate occur until steady state condition again reached asymptotically. 
shown Fig. this phenomena satisfactorily simulated both the 
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dynamic wave and proposed diffusion wave models but ignored completely 
the kinematic wave model. 

The computations for the two examples presented here was performed 
IBM 370/145 computer the Middle East Technical University. For each 
one the examples, the execution time was approx sec, 150 sec, and 480 
sec, respectively, for the diffusion, kinematic, and dynamic wave models. These 
values are presented merely reference. The noticeably different execution 
times for the kinematic and the diffusion wave models are partly due 
smaller computation grid adopted the former, and partly due faster 


0 
Time, min 


FIG. 9.—Comparison Results Example (a) Channel Channels and 
Channel and (d) Channel 


convergence numerical results the latter model. The difference computer 
execution times for these two models smaller for other examples (not reported 
here) having identical grid sizes and numerical schemes which show consistently 
shorter execution time for the proposed diffusion wave model than the nonlinear 
kinematic wave model. The difference the execution times for the diffusion 
and the dynamic wave models arises from faster convergence the former 
and the difference the structures the two models and expected increase 
exponentially with increasing size the network considered. 
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nonlinear diffusion wave model for flood routing dentritic-type open-chan- 
nel networks developed and presented this study. Comparison this model 


with dynamic wave and nonlinear kinematic wave models leads the following 
conclusions: 


The proposed diffusion-wave model can satisfactorily simulate the mutual 
backwater effects channels joining junction. 

The proposed model nearly accurate the dynamic wave model 
which may classified among the most sophisticated one-dimensional routing 
techniques available the literature. 

The proposed model faster and cheaper computation than kinematic 


wave model which may considered one the simplest hydraulic routing 
models known. 
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The following symbols are used this paper: 


flow cross-sectional area; 
incoming channel; 
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incoming channel; 

out-flowing channel; 

any reach function; 

function flow parameters; 

residual function 

gravitational acceleration; 

head loss through junction; 

water surface elevation; 

initial guess for 

improved guess for 

index for flow sections; 

number flow sections channel dynamic- and kinematic-wave 
models; 

index for time stage; 

constant; 

index for incoming channels; 

number flow cross sections along channel 
number flow cross sections along channel 
number flow cross sections along channel 
Manning’s roughness factor; 

index for out-flowing channels; 

discharge; 

hydraulic radius; 

channel bottom slope; 

friction slope; 

storage junction; 

time; 

cross-sectional averaged flow velocity; 
horizontal distance; 

correction 

time increment; 

reach length; and 

momentum correction coefficient. 
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STOCHASTIC MODELS SUSPENDED-SEDIMENT 
DISPERSION 


INTRODUCTION 


Considerable progress has been made recent years the analysis dispersion 
dissolved substances turbulent streams (68), and dispersion discrete 
sediment particles moving slowly bed load rolling sliding along the 
stream bed (71). 

Between these two transport modes the class transport typified the 
movement discrete sediment particles carried mainly suspension. Suspended 
particles follow the movement the dispersive fluid, except that they tend 
lag behind and settle under the influence their own inertia and the gravitational 
pull. sizable part the sediment transported natural streams frequently 
carried suspension and when suspension moves rate quite close 
the stream velocity. accurate prediction this suspended load utmost 
importance many instances. This paper solely concerned with the movement 
sediment particles while carried suspension. 

Sediment particles being transported streamflow are suspended turbulent 
fluctuations and further dispersed mean-velocity gradients. addition, the 
mean flow exhibits random characteristics due the irregularities the stream 
boundaries and the daily variations the stream discharge. This behavior 
reflected the variability the suspended load distribution (9). Therefore, 
the transport suspended sediment can justifiably regarded intrinsically 
random process. The conventional theories suspended transport are, however, 
exclusively deterministic (33) and largely restricted the condition steady, 
uniform flow. This limitation prompted number investigators use 
probabilistic methods order derive time and space distributions suspended 
particles under more general conditions. Two distinct approaches developed 
independently over the years. One based the theories turbulent diffusion 
continuous movements, whereas, the other arose from the random walk 
theory stochastic processes. Review these theories’ development and 
applications the main objective this paper. 


Hydr. Engr., U.S. Dept. Agric. Science and Education Administration, 
Sedimentation Lab., P.O. Box 1157, Oxford, Miss. 38655. 
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Evolution Probabilistic theory diffusion continuous 
movements was originated Taylor (79), who considered the Lagrangian history 
single fluid particle homogeneous isotropic turbulent flow. Subsequently, 
Richardson (65), Batchelor (4,5,6), Taylor (80), and others (29) evolved the 
original ideas into more complete theory. Tchen (81) used this approach 
study the diffusion rigid particle carried turbulent flow. doing 
extended the equation motion submerged rigid particles developed 
Basset (3), Boussinesq (10), and Oseen (58) (BBO equation). Corrsin and 
Lumley (22), among others, further developed Tchen’s formulation. Because 
the nonlinear nature this equation, solution not readily obtainable. 
This led many linearizing assumptions subsequent efforts solve the 
BBO equation (35,36,46,72,81). All these studies have introduced strong simpli- 
fications that rendered the solutions inapplicable. other words, all these 
solutions suggest diffusivity the suspended particles equal the diffusivity 
the surrounding fluid, which experimentally incorrect (41,75). Subsequent 
analyses Friedlander (30), Csanady (23), and Peskin (61) introduced more 
realistic approximations and were able show that the difference diffusivities 
stems from the fact that, general, the Lagrangian correlation turbulent 
fluid velocities and the correlation velocities experienced the particles 
not coincide. Peskin (61) went derive conditions under which the 
fluid and particle diffusivities may taken the same. This criterion 
significant the experimental determination sediment diffusion from 
measurements turbulence characteristics (57). 

The random walk concept was introduced Pearson (60), although Rayleigh 
earlier considered formally identical problem while studying the composition 
vibrations unit amplitude with phases distributed random (64). This 
concept became the object many investigations, which were later reviewed 
Chandrasekhar (14) and Kac (39) among others. The random walk theory 
has been applied the solution diffusion problems two distinct manners. 
The classical approach consisted obtaining the solution complex differential 
equations the limiting solution equivalent discrete probabilistic process. 
The classical example the solution the Fokker-Planck diffusion equation, 
which formally equivalent the Kolmogorov the theory 
random functions (31). The other approach direct simulation the 
physical process imitating the behavior representative discrete particles 
through computer simulation. The effectiveness this method was demonstrated 
Bugliarello and Jackson (12,13), Sullivan (77), Chiu (16), Bayazit (8), Sumer 
(78), and and Shen (48). 

This paper was written the instance the ASCE Task Committee 
the State the Art Stochastic Hydraulics. intended review 
some the more significant findings that evolved from the above analyses, 
and presents some suggestions regarding their extension complex problems 
practical interest. The underlying theories are only briefly outlined. detailed 
presentation can found the treatise Monin and Yaglom (56). 


Continuous Movements 


Diffusion Stationary Homogeneous Turbulent Flows.—Consider set 
discrete rigid particles diffusing statistically steady homogeneous turbulent 
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field without mean shear. The particles are all supposed have the same settling 
velocity and fill the volume time and all subsequent times 
(Fig. 1). Assuming that the particles are initially distributed with uniform 
concentration C,, the problem find their concentration distribution 
subsequent times. This problem equivalent determining the probability 
that particle lies within the space volume given that 
the initial probability the particle being within (4). 

Assume that each particle not influenced the neighboring particles, that 
is, the distribution its displacement independent the other particle 


FIG. Suspended Particle 


trajectories. the conditional probability particle displacement 
given its initial position x,, then 


which the integration over all space. This equation can also written 
the convolution 


The Fourier transform this convolution yields the following characteristic 
function 


Assume further that the probability distribution the particle displacements 
Gaussian with nonzero mean given w(t The characteristic function 
for the displacement then 
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x,, and [-] denotes average over the ensemble all particles. Eq. 
substituted into Eq. and, after computing the inverse transform, the result 
differentiated with respect time. This yields 


P(x, 


Noting the relation the probability the concentration point 
Eq. can rewritten 


2 


which This derivation based the work Batchelor (3), 
which was extended Barfield, al. (2) the case external gravitational 
field. 

Consider now homogeneous turbulent flow with uniform mean velocity 
the direction. the absence gravitational field, similar derivation 
assuming Gaussian displacement around mean U(t gives 


Eqs. and are particular cases the equation for diffusion suspended 
sediment that can derived from the conservation mass principle (70). 

The time rate change the displacement covariance matrix appearing 
Eqs. and can written (4,6) 


0 


0 

which the instantaneous Lagrangian particle velocity; and the 
particle Lagrangian velocity correlation tensor. Experience has shown that the 
Lagrangian correlation function isotropic turbulent velocities has shape 
similar the Eulerian correlation function (35). Assuming that this holds true 
for the solid particle correlation and considering the turbulent motion 
isotropic, R,,. addition, for small diffusion times. Thus, 
integrating Eq. twice results the displacement variance 


the diffusion time increases, the integral Eq. approaches asymptotically 
the Lagrangian time integral scale 
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and Eq. yields the limit 


Consequently, the variance particle displacement varies for small diffusion 
times, whereas for diffusion times the variance varies linearly with 
time and the diffusion process follows Fick’s law (i.e., diffusion with constant 
transfer coefficient). The latter behavior direct consequence having assumed 
the particle displacements normally distributed around the convective mean. 

The preceding relationships have been derived describe diffusion 
stationary, homogeneous turbulent field without mean shear. Thus, they would 
seem inapplicable diffusion open channel flows finite cross sections. 
Moreover, these cases the longitudinal differential convection combined with 
lateral turbulent mixing gives rise longitudinal dispersion (27), usually more 
important than longitudinal diffusion. However, Batchelor and Townsend (7) 
pointed out that steady, uniform flows confined rigid boundaries the 
turbulence homogeneous and nondecaying the streamwise direction, and 
the velocity fluid particle stationary random function time about 
the mean flow velocity. Although confined stream nonhomogeneous 
the cross-stream direction, wandering particle moves randomly towards one 
wall and then towards the other and its velocity becomes statistically stationary 
after long enough period time for the influence the particle’s initial 
position vanish. Under these conditions, Eq. applies longitudinal 
dispersion uniform shear flows. Taylor (80) observed that far enough down- 
stream from the source the concentration gradients associated with the velocity 
shear would balanced the tendency the lateral turbulent mixing 
decrease those concentration gradients. found that the spread the cross- 
mean value concentration, relative coordinate system moving 
with the mean flow velocity, produces indeed apparent Fickian diffusion 
along the channel axis. Elder (25) extended Taylor’s asymptotic analysis 
dispersion infinitely wide open channel. More recently Fischer (27) gave 
method for applying Taylor’s analysis natural streams. Also, Fischer (28) 
noted that far Taylor’s theory concerned, stream curvature and 
nonuniformity can tolerated the reach long enough, and provided that 
the nonuniformity random such way that the Lagrangian velocity remains 
stationary. The preceding observations support the applicability 
Fickian diffusion theory open channels for large dispersion times. This 
conclusion has been confirmed number experimental investigations 
(26,57,70). They all restricted their measurements the principal components 
the diffusivity tensor which not depend the assumption local isotropy. 
Some the results these investigations are summarized below. 

Sayre and Chang (70) measured the longitudinal dispersion suspended 
silt particles wide open-channel flow having average velocity 
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0.35 m/s and depth 0.24 Concentration-versus-time distributions 
were measured several dispersion distances collecting samples rapid 
succession the dispersant passed sampling station. The distributions measured 
mid-depth not differ significantly from the distributions measured 1/4 
and 3/4 the flow depth because appreciable settling occurred over the 
range dispersion distances. the diffusion time increases the concentration 
distributions approach Gaussian shape, suggesting that the fine silt diffuses 
accordance the normality hypothesis. The time variances these distribu- 
tions, are plotted against the travel distance, x,, Fig. Although there 
considerable scatter, the data tend follow the linear relationship between 
and derived Yotsukura (Ref. 83; see also Eq. 11, Ref. 70) for 
one-dimensional Fickian diffusion turbulent open channel flows: 


Here D,) constant diffusion coefficient. The value used 
Fig. was determined Sayre and Chang (70) from dye dispersion experiments 
for the same flow conditions used the silt dispersion studies. 
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FIG. 2.—Time Variance Function Travel Distance (After Ref. 70) 


Engelund (26) conducted dispersion experiments where tracked the motion 
discrete particles wide rectangular open channel flow 0.055 deep 
with average velocity 0.30 m/s. Buoyant spheres diam were 
released the free surface and their longitudinal, and lateral, displace- 
ments were measured subsequent diffusion times. After repeating the experi- 
ments many times, determined the variances those displacements. The 
corresponding rms velocities and time integral scales were obtained from 
set turbulence measurements. The results, shown Fig. corroborate Eq. 
12. The surface turbulence was found anisotropic, the longitudinal diffusion 
coefficient being nearly twice the lateral coefficient. 
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Eq. describes the turbulent diffusion terms Lagrangian characteristics 
particle moving along its path. practice, the Lagrangian characteristics 
are not easily obtainable and most the information available consists Eulerian 
descriptions. Thus, important relate these two types information 
order express the Eulerian results Lagrangian form. homogeneous, 
steady, turbulent flows, the Eulerian time-averaged turbulent velocity, 
commonly assumed equal That is, sampling many different particles 
they pass fixed point taken equal instantaneous ensemble average 
many point velocity traces. Lumley (50) presents formal arguments which 
validate the equivalence these measurements under the restrictions homo- 
geneous, isotropic, and stationary turbulence. 

The problem estimating turbulent diffusivities from point measurements, 
therefore, becomes how approximate the Lagrangian autocorrelation function 


e@eece Diffusion (Ref. 26) 


FIG. Variance Function Travel Time (After Ref. 26) 


and its integral scale. Hay and Pasquill (34) proposed that the Eulerian autocor- 
relation might transformed into the Lagrangian autocorrelation. That is: 


(t) 


Here dimensionless constant related the time integral scale ratio 


which the Eulerian time integral scale. Although subsequent work tended 
support this approximation (73), became apparent that was dependent 
other turbulent characteristics (1), invalidating the empirical constant value, 
proposed Hay and Pasquill (34). Because the convenience afforded 
Eq. number studies sought clarify the functional properties 
most notably that Philip (62). This analysis was based the 
hypothesis proposed Corrsin (19,20). The Lagrangian autocorrelation function 
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formed from ensemble average the time-lagged product single-particle 
velocities having common origin. According Corrsin’s hypothesis this function 
can estimated from the general Eulerian space-time correlation function relative 
coordinate system moving with the mean flow velocity, properly sampling 
(or weighting) the Eulerian data account for the particle distribution space 
and time. Then 


interval (x,x dx) time and the origin coincident with the mean 
convective frame reference. 

Philip (62) assumed homogeneous, isotropic, stationary turbulent motion, 
and introduced the functions 


2 
x 


the Eulerian space-time correlation integral scale 


9 


0 0 


which results from integration Eq. Combining Eqs. 16, 17, 18, and 20, 
Philip obtained integral equation for that solved iteration for 
various values the Eulerian diffusivity parameter 
Introducing the results Eq. yielded, after numerical integration, the 
Lagrangian-to-Eulerian space-time scale ratio 

which turned out decaying function a,, that borne out 
data reported Baldwin and Mickelsen (1). The behavior this function 
such that the Eulerian parameter approaches zero, while 
F(a,) becomes approximately inversely proportional for 


express the results terms Philip (62) employed Taylor’s hypothesis 
over the entire time span, namely: 
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u 


supported available atmospheric and laboratory data (see Fig. Ref. 
62). Clearly, Eq. does not predict unique dependence turbulent 
intensity. 

Knowledge the time integral scale ratio permits the evaluation the 
Fickian diffusion rate from Eq. using Eulerian measurements. Nordin and 
McQuivey (57) used this approach compute the vertical distribution mean 
sediment concentration two-dimensional steady uniform flow. this case, 
Egs. 12, and yield 


Integrating yields 


which the concentration some reference level The values 

and measured McQuivey, al. (53) were used evaluate 

the transfer coefficient. The results compare favorably with the values obtained 


Jobson (37) from dye diffusion measurements. Jobson’s data follow closely 
the parabolic distribution 


constant. can shown that this constant takes the value two-dimensional 
turbulent flows, where the von Karman’s constant and the bed shear 


which w/xu,. This equation was first introduced Rouse (66) and 
has been the base for many deterministic studies suspended sediment load 
(33). 

The preceding example illustrates the combination probabilistic and deter- 
ministic methods obtain practical solutions. The probabilistic theories permit 
point estimates local diffusion rates, while solutions the mass transfer 
equation incorporate the appropriate boundary conditions well the sediment 
settling. Other illustrations this combined approach are discussed Sayre 
(67,69). 

Relative Diffusivity Solid and Fluid Particles.—The foregoing results assume 
that the solid particles follow exactly the motion the surrounding fluid. 
shown the following, this leads the experimentally incorrect conclusion 
that the particle diffusivity identical, general, the fluid diffusivity. 
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thus important investigate the conditions under which the foregoing 
assumption valid. 

The problem relative diffusivity was first considered Tchen (81). His 
analysis, and subsequent clarifications others (22,35) extended the BBO 
equation the case spherical particle moving turbulent flow. The 
particle Reynolds number taken small enough that the particle drag 
can described Stokes’ law. The particle regarded part sparse 
cloud diffusing turbulent domain infinite extent, order neglect 
the effect boundaries and particle interaction. After assuming that the particle 
encounters uniform’’ flow along its trajectory, and neglecting 


all nonlinear resistance, Tchen (81) obtained the following first-order linear 
equation 


the solid particle velocity along its trajectory; u,, the ith velocity component 
the uniform’’ flow; the solid particle diameter; the specific 
gravity the particle relative the fluid; the fluid kinematic viscosity; 
the time the beginning the particle trajectory; the ith component 
the gravitational field force. Tchen (81) used Eq. obtain the particle 
diffusivity terms the stationary Lagrangian energy-spectrum 


function. Chao (15) further generalized Tchen’s analysis and related studies 
(23,30,74) treating Eq. stochastic equation subject Fourier 
transformation. dropped the deterministic external force and introduced 
expressions defining the Fourier transform the fluid and solid particle velocities 
and the corresponding energy-spectrum density functions, and Using 
the Wiener-Khintchine theorem, the stationary autocorrelation coefficients for 
the fluid and solid particle velocities were expressed Fourier cosine transforms 
the energy-spectrum densities. Following Tchen (81), Chao (15) assumed 
that small solid particles would tend follow the fluid motion and that the 
Lagrangian correlation would thus serve good approximation the particle 
correlation function. This also the case for which the assumption Stokesian 
motion the particle most likely hold. Introducing the Lagrangian correlation 
functions Eq. Chao (15) formulated the ratio the solid particle, 


p? 
fluid, D,, diffusivities 


0 


which defines the harmonic frequency component. This ratio was shown 
Chao approach unity large diffusion times, regardless the particle 
inertial properties and the structural features the flow. This result, also obtained 
Tchen (81), not supported available experimental evidence (41,75). 


SUSPENDED-SEDIMENT DISPERSION 743 


This disagreement stems from the assumption that the solid particle always 
the same fluid elements its neighborhood. However, the particle 
made follow the irregular path the fluid accelerations induced drag 
forces exerted the fluid itself. These drag forces can only exist the expense 
relative velocity between the particle and the fluid. One effect this 
relative motion that the particle forced continuously change its fluid 
neighborhood and therefore lose its velocity correlation more rapidly than 
the fluid. This has been termed the crossing Yudine 
(84). Consequently, better definition solid particle diffusivity should contain 
some measure the time-space correlation between the particle and fluid 
velocities. significant step this direction was taken Peskin (61). 
Considering heavy particle, Eq. may approximated 


and the ith component the Eulerian fluid velocity, encountered 


Solid Particle Path 


Fluid Element Path 


FIG. 4.—Effect Crossing Trajectories 


the particle. The turbulent velocity will considered random function 
with zero mean relative some uniform mean flow velocity. Integrating Eq. 
gives (49) 


which the initial particle position has been omitted for convenience. This 
nonlinear stochastic equation since the velocity random function 
the random particle displacement. Peskin (61) solved indirectly statistical 
sense, considering the most probable fluid velocity encountered the particle 
given position. 

Fig. shows the trajectories solid and fluid particles coinciding 
Because particle-fluid slip occurrence, the particles’ Lagrangian displace- 
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ments and are not necessarily the same later time While the fluid 
element has Lagrangian velocity the fluid velocity encountered 
the solid particle will one many possible velocities associated 
with fluid particles having positions ..., Estimation particle 
diffusivity requires relating the velocity encountered the solid time 
the velocity the fluid element coinciding with the solid particle 
This relation obtained computing the most probable fluid velocity 
encountered the particle given certain fixed Lagrangian fluid process. 
Mathematically, this probable velocity expressed the conditional expectation 
are distributed according joint Gaussian distribution with Eulerian correlation 
gives 


The equation motion (Eq. 31) can now written 

0 


This equation replaces the true trajectory the particle the path 
average particle depending stochastically the given fluid Lagrangian informa- 
tion. further assuming the turbulence isotropic, and restricting 


the relative displacement small values, the correlation can 
expanded 


which the Eulerian length microscale the turbulent field. Using Eq. 


33, expression for obtained, which squared and substituted back into 
yields 


1 s s 


and The mean square value the particle displacement 
obtained integration over the joint probability density the Lagrangian 
fluid velocities v(s), etc. This results integral equation requiring 
information the fourth-order fluid-velocity moments. closure assumption 
introduced taking quasinormal form for the fourth-order moment, and 
using the following form the velocity Lagrangian correlation (42) 


L 
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The final expression for the mean square displacement large times (61) 
and from 


2 


Introducing now the Lagrangian length integral scale (79), and 
dividing Eq. the turbulent fluid diffusivity, yields the following 


diffusivity ratio: 
3K? 
K+2 


Eq. plotted Fig. along with data measured Kennedy (41) and 
later corrected Peskin (61). Both the inertial parameter and the turbulence 
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FIG. 5.—Rate Particle Diffusivity Fluid Diffusivity (After Ref. 61) 


scales have significant influence the relative diffusivity solid particles. 
physical interpretation this relationship not difficult. When flow 
highly correlated the fluid particles the neighborhood 
the solid particle all have similar velocities. Consequently, the correlation 
fluid velocities encountered the particle will similar the Lagrangian 
fluid correlation and the diffusivity ratio will approach unity. the contrary, 
becomes larger than the particle will almost random fluid velocities 
and its diffusivity will decrease relative the fluid. 

most turbulent flows the length scale comparable the lateral Eulerian 
integral scale This scale associated with large-scale eddies; thus its magnitude 
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may taken equal the flow depth two-dimensional flows. The ratio 


microscale always much larger than the integral scale Therefore, one 
must have for the diffusivity ratio close unity. Introducing 
and 0.02U, Eq. yields for sand particles diffusing water 


which AU/v the flow Reynolds number. obvious that this 


condition satisfied the experiments Engelund (26) and Sayre and 
Chang (70). 


Consider discrete particle undergoing sequence random independent 
displacements magnitude Ax, the direction. Let denote the probability 
one-step transition the positive direction, and the probability one-step 
transition the opposite direction. Further assume that the process stationary 
and Markovian. Then the probability the particle being time after 
displacements satisfies the difference equation 


Expanding means Taylor’s theorem gives 


neglecting terms order and When passing the limit, 
appropriate ratio Ax, must maintained the process not degenerate 
During time the total displacement the sum independent steps having 
mean q)Ax, and variance these values are remain 
These conditions are satisfied 


which and suitable constants. The mean and variance the net 
displacement time are now given 2ct and respectively. and 
are interpreted the mean convective velocity, and the turbulent diffusion 
coefficient, respectively, Eq. becomes 


which known the Fokker-Planck equation. This result has been generalized 
Kolmogorow multidimensional equation with variable coefficients (31). 
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the simulation applied the dispersion large number particles 
the probability multiplied the total number realizations gives the 

Eq. shows that random walk simulation replaces, the limit, the solution 
the Fickian diffusion equation. Taylor (79) recognized this introducing 
the Lagrangian integral time scale. This integral can regarded measure 
the interval over which the particle velocity retains correlation with its 
initial value. the duration the random walk steps are selected such that 
T,, then the particle displacement time may represented 
sum small displacements that are weakly dependent random variables. 
Therefore, according the central limit theorem, the probability distribution 
the sum will approach the normal form the number steps becomes 
very large. Furthermore, Goldstein (32) has shown that random walk model 
embodying the behavior turbulent diffusion small and large times leads 
the so-called For small times this equation reduces 
the wave-propagation equation; for large times becomes diffusive equation. 
Corrsin (21) has shown that this dual behavior analogous the asymptotic 
behavior turbulent diffusion (Eqs. and 11). These observations serve 
illustrate the fact that the random walk and the continuous movements theories 
are basically equivalent. 

Examples some typical applications the random walk technique its 
more basic form, direct simulation, the study suspended-sediment dispersion 
follow. 

Settlement Solid Particles Two-Dimensional Flows.—Bayazit (8) developed 
model akin random walk process simulate the motion solid particle 
settling steady uniform open-channel flow over smooth bed. Incorporated 
the model are path lengths, fluctuating velocity scales, mean velocities, etc., 
that are functions position. The particle assumed follow exactly the 
motion the surrounding fluid moves downwards through the statistical 
field variables with settling velocity Once the particle reaches the stream 
bed, the motion terminated. This occurrence defines the settling length 
the particle. Many such realizations are used compile the probability 
distribution the settling lengths. 

Because the lack available Lagrangian information, attempt made 
simulate the process using the available Eulerian information pseudo-La- 
grangian way. The turbulent velocity fluctuations the longitudinal and lateral 
directions are neglected, and the vertical fluctuations are assumed normally 
distributed around zero mean. That 


2 

which the rms obtained from experimental data analyzed Partheniades 
and Blinco (59). The vertical fluctuation assigned certain value each 
time level means random number generator. During each time step 
the particle experiences horizontal convection displacement Ax, 
followed randomly directed vertical displacement Ax, (u, 
where the local mean velocity. The time step selected render 
equal the turbulence Eulerian microscale estimated from the data collected 
Raichlen (63) and McQuivey and Richardson (54). 
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This model was later extended and Shen (47), who introduced several 
refinements. They computed the local mean velocity from the velocity-defect 
law, which valid over larger depth range than the law the wall used 
Bayazit (8). They further assumed that the streamwise fluctuations are not 
negligible and normally distributed, such that 


The rms computed from 


2 


which the covariance obtained from simplified form the Reynolds 
equations governing the two-dimensional mean flow; correlation coefficient 
derived from available experimental data (43,54); and defines the ratio (54) 


1/2 


the light Eqs. and 48, and Shen assumed and have 
bivariate normal joint distribution. Thus, they obtained 


c 


each time step the rms and are computed from and 50; 
then and are randomly generated using and 51. The time step and 
vertical displacement are computed following Bayazit (8), and the convective 
displacement obtained from At. However, contrary Bayazit, 
not treated random quantity, but local mean convective time 
increment chosen that the expected value |Ax,| equal the local 
Eulerian length scale. 

The random walk model was simulated computer for various flow 
conditions and particle characteristics. each case the random walk was realized 
more than 100 times, and the mean value and standard deviation the settling 
length distribution were computed. The computed values compare favorably 
with the measurements Bayazit (8) and Lean (45), shown Fig. The 
frequency histograms the simulated settling length distributions for sediment 
mixtures agree with the experimental results obtained Jobson and Sayre 

the preceding models the solid particles move between reflecting boundary 
the free surface, and absorbing boundary the channel bed. similar 
model using reflecting barriers both the bed and free surface was developed 
Sumer (78). This model was basically patterned after those Sullivan (77) 
and Bayazit (8). 

Dispersion Sediment Three-Dimensional Turbulent (16) was 
the first formulate random walk simulation solid particles dispersing 
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three-dimensional turbulent flow. His analysis leads three-dimensional 
analog Eq. 44. the preceding examples, the random walk step considered 
the resultant convective component (which incorporates the particle 
settling), and random diffusive component. Because the lack information 
the distribution turbulence characteristics three-dimensional flows, Chiu 
represented the diffusive step random vector. The model also assumes 
that there particle-fluid slip, and that interactions between particles 
not exist. The time increment determined equating the diffusion coefficient 
momentum transfer coefficient based von Karman mixing length model. 

Chiu and Chen (17) applied this model simulate the dispersion 


Bayazit (8) 
Lean (45) } Mean Settling Length 


Bayazit >) Standard Deviation 
Lean (45) J of Settling Length 


= 
2 
> 
2 


Measured Values (cm) 


FIG. Simulated and Measured Settling Lengths (After Ref. 48) 


sand particles steady, uniform open channel flow. Measurements conducted 
0.60-m deep flow with average velocity 1.10 m/s yielded information 
the distribution longitudinal mean velocity and sediment concentration 
cross section located the established flow region (18). The sediment 
distribution over the cross section was computed continuously releasing sand 
particles from hypothetical point source. This source was located upstream 
near lower corner the channel, where the presence the side wall and 
the corner eddies render the mean flow three-dimensional. The transversal mean 
velocity components were computed from the measured velocity distribution 
using the method developed Ligget, al. (48). The simulation consisted 
releasing large number particles, one time, and tracking their 
trajectories. This equivalent releasing all the particles simultaneously because 
the stationarity the flow. some distance from the source the influence 
the source location will cease felt the particles. Thus, the concentration 
will approach the same steady distribution independent whether the particles 
are being released from point distributed source. Fig. presents the 
simulated steady state concentration distribution, the form lines equal 
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concentration, along with the measured values. this figure the concentration 
values are normalized the total amount sediment the entire cross section. 
Although the data are quite limited, they compare well with the simulated values. 
The simulation reflects the effects gravity, resulting increase 
concentration towards the channel bed. 

The preceding experimental verifications tend substantiate the random walk 
simulations reviewed this paper. However, there need for additional 
refinements make this technique adaptable more general situations. Some 
aspects susceptible further development are now considered. 

The parabolic character Eq. implies that sediment particles introduced 
fluid stream will spread instantaneously throughout the length the stream. 
This clearly contradicts the finiteness velocities associated with actual fluid 
flows. related the fact that only the concept mean, and not instantaneous, 
particle velocity valid when only the particle displacement interpreted 
Markovian function (56). This may lead some cases unacceptable errors 


Simulated Concentration 


; Measured Distribution 


Values 


FIG. Simulated and Measured Steady-State Distribution Sus- 
pended Sediment (After Ref. 17) 


close the edges the dispersing cloud. Therefore, desirable have 
models based more general theory which the boundedness the particle 
velocity would taken into account. Monin (55) proposed one such theory 
were not the particle displacement but the two-dimensional random function 
formed the particle position and velocity Markovian. 

There also need for more general dispersion models capable handling 
strongly nonstationary flow fields. They will find application instances where 
the interest focussed the dispersion sediment relatively short channel 
reaches irregular geometry. These models should free the approximations 
made the existing models regarding the turbulence structure, and should take 
into consideration the dependence solid-particle diffusivity the turbulent 
flow features well particle inertia. promising approach was presented 
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Kau and Peskin (40). They used improved version Deardorff’s 
Navier-Stokes solver (24) investigate the Eulerian structure three-dimen- 
sional turbulent flow. Mean flow characteristics and Eulerian statistics were 
obtained from the solution. approximate form Eq. was integrated 
obtain Lagrangian information the motion solid particles released 
the streamflow. Lagrangian statistics describing the particle diffusivities were 
obtained from the simulated particle trajectories. This study another good 


illustration effectual combination probabilistic and deterministic tech- 
niques. 


Summary 


Prediction suspended-sediment dispersion natural streams important 
the evaluation sediment deposition patterns and transport sediment-borne 
pollutants. Starting with the early theory diffusion continuous movements, 
and proceeding more recent random walk models, existing probabilistic methods 
for predicting turbulent dispersion suspended sediment were reviewed. 

All these models assume that the sediment suspended concentrations 
that are not large enough significantly affect either the fluid velocity 
the motion the suspended particles. These are all kinematic models the 
sense that they require information the dynamics the dispersive medium. 

The theory diffusion continuous movement can extended natural 
channels provided the Lagrangian mean velocity remains statistically stationary. 
those cases which the particle properties and turbulence characteristics 
result small value the Peskin’s parameter (Eq. 40), the sediment diffusivity 
may obtained from either Lagrangian Eulerian information about the 
turbulence characteristics. 

Markovian random walk models are equivalent Taylor’s theory the limit 
large dispersion times. They are not restricted principle stationary flows, 
but impose probability-distribution restrictions which are not present the 
continuous movements technique. The random walk simulations depend for 
accuracy knowledge the turbulent flow characteristics. These data are 
very limited and mostly restricted two-dimensional flows. This obstacle may 
surmounted, however, using Navier-Stokes solvers when experimental 
information available. These solvers have already been introduced turbulent 
boundary-layer computations, and workers this field are extending them 
include complex three-dimensional flows (11,44,76). The dispersion models 
considered this paper assume that discrete particles follow exactly the turbulent 
fluctuations. This assumption usually limits the permissible particle size range 
silts and fine sands. Nevertheless, this constraint can relaxed coupling 
the flow models suitable form the BBO equation. 

The ultimate evolution probabilistic methods into formulations general 
applicability may still lie ahead. There still large gap between the models 
that have been developed far and the practical solution complex hydraulic 
problems associated with the transport suspended sediments. For the present, 
the combination probabilistic and deterministic solutions seems the 
most efficient approach practical applications, where the strengths both 
methods are used while their individual weaknesses are complemented. Nev- 
ertheless, large-scale simulation sediment dispersion using probabilistic-de- 
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terministic approach cannot regarded practical even with today’s largest 
digital computers. However, such approach can used reduced scale 
complement the fully deterministic models derived from continuous theory. 
The weakest part these deterministic models the lack accurate procedures 
for predicting the effect channel curvature and irregularities along short reaches, 
the spatial variation mass transfer coefficients. This information could 
obtained simulating the dispersion suspended particles fashion 


similar those propounded Chiu and Chen (17), and Kau and Peskin 
(40). 
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ratio rms transverse velocity fluctuations 


streamwise velocity fluctuations; 
turbulent fluid diffusivity; 

solid particle diffusivity; 

solid particle size; 

energy spectrum density functions; 
probability density function; 
scalar component gravitational field force; 
flow depth; 

Eulerian intensity turbulence; 

Peskin’s parameter; 

Fickian diffusion coefficient; 

lateral Eulerian integral scale; 

Eulerian length integral scale; 

Lagrangian length integral scale; 

probabilities; 

flow Reynolds number; 

Eulerian velocity correlation; 

macroeddy Reynolds number; 

Lagrangian velocity correlation; 

specific gravity solid particles; 

Eulerian time integral scale; 

Lagrangian time integral scale; 

Eulerian space-time correlation integral scale; 
time; 

cross-sectional average velocity; 

local mean velocity; 

Eulerian velocity fluctuation; 

mean-square Eulerian velocity fluctuation; 
scalar component Eulerian fluid velocity; 
scalar component Eulerian particle velocity; 
bed shear velocity; 

volume; 

Lagrangian fluid velocity; 

scalar component Lagrangian particle velocity; 
solid particle settling velocity; 

relative displacement; 

vector positions; 

coefficients the BBO equation; 

Eulerian parameter; 

probability distribution particle displacements; 
fluid-solid particle relative displacement; 

time increment; 

displacement increment; 

time scale ratio; 

mass transfer coefficient; 

relative depth; 

von Karman turbulence coefficient; 
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Eulerian length microscale; 

fluid kinematic viscosity; 

wave number vector; 

correlation coefficient; 

constant derived from Stokes drag law; 

time variance; 

displacement variance; 

harmonic frequency component; and 

denotes average over the ensemble all marked particles. 
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RUN FOR PREDICTION 


INTRODUCTION 


For most practical applications, flood sufficiently defined the complete 
flood hydrograph and its exceedance probability known. Important elements 
the flood hydrograph are its instantaneous peak discharge rate, its total flood 
volume, and its shape. The flood peak and its exceedance probability give enough 
information for many practical problems, e.g., the sizing storm drainage 
pipes, some culvert sizes etc., where storage unimportant. There are many 
methods which can used estimate the peaks and their exceedance probabil- 
ities. All these methods are explicitly tacitly based the marginal probability 
distribution flood peaks. Flood volumes can treated similar way. 
Flood volumes are important problems where storage plays important 
part, e.g., the sizing spillways for dams, the calculation flood damage, 
and flood routing. Flood volume analysis very seldom done isolation 
because actual data flood volumes scarce and relatively difficult process. 
Flood volumes alone also have limited value because they not supply 
information depth duration inundation. Therefore, the complete flood 
hydrograph needed. The only general method which widely used, and 
purports supply the complete hydrograph, the unit hydrograph method 
and its derivatives. The unit hydrograph technique, however, suffers two 
important limitations. Firstly, many observed floods are composite floods and 
are therefore not easily described unit hydrograph methods. is, course, 
possible combine design hydrographs, derived the unit hydrograph method, 
into one composite flood hydrograph, but this illustrates the second limitation—the 
return period this composite hydrograph unknown! Indeed, any design 
hydrograph derived through standard unit hydrograph technique, concep- 
tual modeling, has unknown return period. This comes from the combination 
peak and volume whose cross correlation not unity. This serious 
limitation because most analyses floods and designs involving floods require 
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that the return periods the floods known. The importance relating 
flood peaks and flood volumes was also expressly stated Ref. 10. The few 
efforts towards extracting the composite flood hydrographs with the desired 
return period (5,7), were not very satisfactory, which highlights the limitations 
the unit hydrograph technique. 

The problem addressed this paper, therefore, involves two related variables, 
namely flood peaks and their associated volumes. immediately apparent 
that more than one combination flood peak and flood volume may have 
exactly the same probability being jointly exceeded. For example, relatively 
large flood peak combined with small volume may have the same return 
period equally large volume combined with small peak. (Here, and 
what follows, means the reciprocal the probability 
joint exceedance peak and associated volume.) This means that for each 
catchment there exists possible family hydrographs for any selected return 
period. The design hydrograph obtained unit hydrograph technique may 
only one member one the families. 

method which demonstrates the extraction such families hydrographs 
from the flow record river was described previous paper (7). this 
paper, the authors assumed fixed shape for all the hydrographs and they 
described the return periods for the hydrographs terms bivariate lognormal 
probability distribution peaks and surplus run lengths selected base value. 
The volume flood runoff, although not directly modeled, follows from the 
assumed shape the hydrograph. This method was not entirely satisfactory 
because the assumed hydrograph shape looked unnatural and the method resulted 
general underprediction flood volumes. The direct identification the 
flood volume one variable the bivariate description, instead indirectly, 
should improve the method. 

This paper reports the results obtained during the analysis the instanta- 
neous continuous streamflow records selected South African catchments. 
Altogether 1,172 individual flood hydrographs were analyzed, from catchment 
sizes ranging from mile-41,730 mile (26 


Basis 


Bivariate Lognormal Probability Distribution.—The bivariate lognormal proba- 
bility distribution function was selected possible suitable descriptor for 
the flood peaks and their associated volumes. This function relatively easy 
manipulate, and experience has shown that the univariate lognormal often 
satisfactorily describes the probability distribution flood peaks. Use the 
standardized form the bivariate normal function the natural logarithms 
the flood peaks and volumes makes comparisons across many catchments 
possible. This function standardized form (7) 


the two variables and respectively. this application, represents 
the natural logarithm the flood peaks; represents the natural logarithm 
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the associated flood volumes; and represent the means; and 
represents the standard deviations the variables and respectively; and 
the cross correlation coefficients for the variables. 
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FIG. Standardized Bivariate Normal Function with Observed Floods 
Shown and Cross Correlation 0.85 
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FIG. 2.—Bivariate Return Period Contours Positive Quadrant Standardized cdf 
with Cross Correlation 0.85 


Fig. shows the bivariate normal density function for cross correlation 
coefficient 0.85 the form equal probability contour lines. this figure, 
the standardized flood peaks, and their associated standardized volumes, 
for the selected catchments are shown. (Note that data several catchments 
with different are shown; thus the return period contours are only valid for 
those catchments with 0.85.) 
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For practical applications useful present results terms return 
periods. Lines equal return periods can obtained through volume integration 
under relevant parts the probability density function. integration method 
described Tihansky (12) was used. This method seeks points which define 
stipulated volumes under the probability density function within the 
north-easterly space formed when vertical planes are dropped through this point, 
with one plane parallel the y-axis and the other plane parallel the x-axis. 
(An easy way visualize this operation consider the problem dividing 
part cake into equal portions with short knife which does not reach 
the center the cake. Cuts must made parallel the rectangular coordinate 
axes, and the cake circular and symmetrical, the leftover part must again 
circular.) The locus which connects these points then forms contour 
equal return period. Fig. shows such equal return period lines the positive 


TABLE 1.—Bivariate Return Period Contours Positive Quadrant Standardized 
with Cross Correlation Coefficient 0.85 (See Fig. 


Standardized Standardized 
Natural Natural 


Return Return 
Logarithm 


eriod eriod 
contour contour 


(2) (5) 


quadrant the and y-axes system for cross correlation coefficient equal 
0.85. Only the positive quadrant shown because engineers are mostly 
interested floods with both peaks and volumes larger than the mean the 
peaks and volumes respectively. (Table gives the actual values the line 
crossings with the contours.) 

the lognormal bivariate probability distribution function found suitable 
for the description the joint occurrences flood peaks and volumes, then 
Fig. and similar figures for different cross correlation coefficients can 
used for practical applications, demonstrated later this paper. The chi-square 
test was used criterion for the goodness fit. The chi-square test 
sensitive the number class intervals selected, but guide the 


fs: 
Line Volume Peak 
1.040 1.040 1.790 1.790 
0.880 1.142 1.500 1.954 
0.700 1.210 1.170 2.024 
0.520 1.246 0.850 2.046 
20 1.390 1.390 100 2.050 2.050 
1.180 1.522 1.720 2.236 
0.930 1.600 1.340 2.304 
0.680 1.627 0.970 2.322 
1.570 1.570 200 2.300 2.300 
1.330 1.720 1.920 2.494 
1.040 1.795 1.480 2.560 
0.760 1.822 1.070 2.573 
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which represents the number events under consideration. The chi-square 
value given 


i=! 
which represents the observed number events class interval and 
represents the expected number events class interval The number 
the number parameters estimated from the sample. 

The magnitude the base value (the level used decide which floods should 
admitted the sample) was found have strong influence the way 
the data fitted the bivariate lognormal distribution measured the chi-square 
test. trial and error, suitable base value can usually found which ensures 
reasonable fit. Fig. shows this trial and error procedure for the Vaal River 
Standerton, with catchment area 3,187 miles (8,254 The test 
was applied the 95% significance level. The figure shows the density functions 


peak 


FIG. 3.—Trial and Error Solution for Suitable Truncation Level 


for the natural logarithms the flood peaks for various base values; Table 
gives information about the chi-square values obtained for each trial and 
the percentage truncation. (The estimation the parameters the density 
functions from the truncated samples explained later paragraph.) The 
base value which gives the lowest chi-square value for the bivariate lognormal 
function was usually selected, provided was above the return period 
peak, gave reasonable chi-square fit for the peaks and was not unrealistically 
high. For the example shown Fig. trial number was selected which 
truncates 44.7% the probability density function. (The return period 
for peaks was selected lower truncation limit ensure independence between 
successive hydrographs the streamflow record.) 

Out the selected catchments, only one catchment could not satisfy 
the chi-square test for peaks, and only stations could not satisfy the chi-square 
test for the bivariate lognormal distribution. The flood volumes did not fit 
well, but the overall fit was considered satisfactory. The cross correlation 
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coefficients for the logarithms peaks and their associated volumes for each 
the selected stations varied from value with mean 0.78 
and standard deviation 0.12. 

Base Value.—It desirable use much information the streamflow 
record possible. This means using partial duration series flood events 


TABLE Procedure Find Best Truncation Level (See Fig. 


for for 
Trial natural logarithm natural logarithm Bivariate Truncation, 
number peaks volumes percentage 


(1) (2) (3) (4) (5) 


1,698.48 88.2 
1,080.60 81.9 
93.98 76.1 
176.39 70.0 
74.79 
38.75 44.7 
73.64 61.5 


(truncation) 


FIG. 4.—Relations between Distributions Total Populations and Truncated Popula- 
tions 


with relatively low base value instead annual series floods. This 
also implies that the return period, years, not simply equal the reciprocal 
the exceedance probability. Using partial duration series floods results 
truncation the population probability density function, with the truncation 
given the level below which the flood events are ignored. Fig. 4(b) shows 
this truncated probability density function. Integration the truncated density 


764 
7.01 
5.91 
6.06 
3.43 
6.59 
4.16 
2.97 
c 
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function results the cumulative distribution function (cdf) shown Fig. 4(a) 
curve SD, whereas the integration the full (untruncated) density function 
results line ABC the cdf. 

The curve the conditional cdf formed normalization. Note that for 
very big floods the cumulative probabilities for both the partial duration series 
(given SE) and the annual series (given ABC) approximately coincide. 

fairly easy correct the partial duration series return periods that 
the return periods are given years. assuming that the arrival flood 


events, which exceed the truncation level, Poisson process, can shown 
(4) that 


which represents the cdf the annual series events; F,(x) represents 
the cdf derived from the associated partial duration series events; and 
represents the mean number events which arrive per year the partial 
duration series. The exceedance probability given the survival function 


and the return period years 


The assumption that the arrival flood events above suitable base value 
Poisson process, seems acceptable for South African catchments 
(5) and fact, further confirmation for this assumption provided Cramer 
and Leadbetter (3). 

For application the bivariate exceedance probabilities consider Fig. Note 
that the shaded area Fig. corresponds the univariate exceedance probability 
for variable Applying Eqs. these exceedance probabilities therefore 
all that necessary relate the bivariate partial duration exceedance probabilities 
the annual exceedance probabilities and thus the return periods, years. 
This becomes more clear later section the paper when the application 
the theory demonstrated actual catchment. 

The truncation the streamflow records also complicates the parameter 
estimation for the bivariate normal probability distribution. The problem is, 
namely, estimate the untruncated population parameters from truncated 
sample. (The writers are the opinion that this complication well worth 
the effort because more realistic fit the large events the selected probability 
distribution function usually results. This means slight emphasis the large 
events truncation the small events.) 

Parameter Estimation.—The necessary parameters for the standardized bi- 
variate normal distribution function are the means and standard deviations 
the natural logarithms the flood peaks and their associated volumes, respec- 
tively, well their cross correlation coefficient. relate the partial duration 
series return periods the annual return periods need, addition, the 
mean number flood events per year well the percentage truncation. 
This brings the number parameters estimated from the data total 
for each catchment. 
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desirable select the base value such way that successive hydrographs 
are independent from each other. Serial independence greatly simplifies inferences 
about probabilities occurrences hydrographs. was found the senior 
writer (5) that consecutive surplus run lengths horizontal truncation levels 
tend become less serially dependent with higher truncation level. base 
value about 1.5-yr return period peak, these surplus run lengths are already 
serially uncorrelated. With this guide, and the desire get reasonable 
fit the data the bivariate lognormal distribution, base values can selected. 
other words, for each catchment base value can selected which ensures 
that the hydrographs which protrude through are almost serially independent 
while the same time giving reasonable good fit the bivariate lognormal 
distribution peaks and volumes. This level called and its logarithm 
(This juggling base values merely satisfy the chi-square test, perhaps 
not justified, due the small size the samples generally available. base 
value can fixed at, e.g., the 1.5-yr return period peak value for all the 
catchments under study.) 

With the base value selected, follows that flood hydrograph can defined 
the locus described the discharge rate between crossing and 
down crossing this base value. This locus then extended from the base 
value the zero flow level. The hydrograph peak the highest point this 
curve. The area under this curve gives the volume the flood. effort 
was made separate the baseflow from the hydrograph, because the baseflow 
forms part the flood and must included. Furthermore, this stage 
was not envisaged predict run hydrograph families from the causative rainfall 
and the inclusion baseflow therefore does not present problem. For the 
calculation the total flood volume, both the rising limb and the recession 
limb the hydrograph had extended from the truncation level zero 
flow level. Various methods were tried and eventually satisfactory method 
was adopted. For this method, weighted average slopes were calculated for 
the rising and recession limb the hydrograph, respectively. Slopes were 
weighted accordance with the lengths portions the limb, above the 
truncation level, which exhibits more less the same slope. The rising limb 
was then extrapolated downwards its weighted slope from the first point 
below the truncation level where reversal slope occurs. The recession limb 
was likewise extrapolated. For most the flood hydrographs the 
base value was relatively low discharge rate, compared its peak, 
that this rather arbitrary extension the hydrograph limbs zero flow level 
could not have large effect the total flood volume. was thus possible 
extract all the peaks and volumes for all the hydrographs exceeding q,, 
objective manner automatic computation. 

The peaks, and the volumes, for the truncated streamflow record 
for each catchment are therefore available. From this data, parameters must 
estimated for the untruncated population. Aitchison and Brown (1) explained 
how the maximum likelihood estimates these parameters can obtained 
with the help specially constructed tables. They also referred iterative 
procedure, based the Newton-Raphson method, which can accomplish the 
same objective without the use tables. 

Assume that there are pairs observations, u,, v,, with ... 
which the sequence logs flood peaks exceeding the 


a 
| 
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log the truncation level; and the logs the associated volumes. The 
estimation procedure summarized follows, following Cohen (2): 


the standard normal density function and 


« 


€ 


the normal survivor function. terms these functions 


From the data calculate 
= (u, u,,)* 


and then estimate from the implicit formula 


which can achieved any convenient root-finding technique such bisection. 
Then calculate directly 


i, = 


giving estimations the population parameters the logs the peaks. However, 
get the joint distribution peaks and volumes the following necessary. 
Calculate the sample statistics 


(14) 
ue (u,— a)’ 
(16) 
-2 (v, ¥)’ 
— S,) 
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The five parameters for the untruncated population which follows the bivariate 
lognormal probability distribution can therefore estimated from truncated 
sample. 

The remaining parameter, namely the average number events which occur 
per year the truncated sample, easily obtained counting. Note that 
found iteration, the sense that the best fit the bivariate lognormal 
distribution, defined the chi-square test, fixes the optimum 

With the flood peaks and their associated flood volumes satisfactorily described 
probability distribution function, only remains describe the time 
distribution runoff rates for each pair. 

The Pearson Type III function, with which engineers are familiar, has shape 
that resembles the average flood hydrograph shape. This function also more 
flexible than two parameter functions and was previously used describe 
hydrographs Yevjevich (13) and Reich (9). For the purpose this paper, 
was decided follow Reich closely. fit hydrograph shapes, the Pearson 
Type III function can expressed 


which represents the flow rate any time represents the peak 
discharge rate; represents the time from the beginning the hydrograph 
the peak; and represents the time from the peak the centroid the 
Fig. shows equation and its parameters. Integration Eq. 
results the flood volume 


i G m/G m 


which the gamma function. According Eq. 26, when for particular 
flood, the volume, and the peak, are known, then and are uniquely 
related. Consequently, only one extra parameter besides and necessary 
describe hydrograph uniquely, namely either The parameter 
this end, various forms linear and nonlinear regressions were performed, 


HY6 
Then, defining 
Sy 
the estimates the remaining population parameters can computed from 
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without much success. Eventually, was decided simply regress 
and and ignore catchment properties. This means that the original 
three-parameter function selected describe the hydrograph shape, degenerates 


two-parameter function, but considerable simplification applications 
results. The final equation was 


with 0.72. 
Here, cubic meters, cubic meters per second, and hours 
35.31 ft). This equation based 1,112 hydrographs with positive 


few hydrographs with negative were excluded from this part the 
analysis. 


For properly-gaged catchments, all the necessary parameters can therefore 
estimated. With and thus, chosen for desired hydrograph, 


centroid 


time 
FIG. 5.—Hydrograph Shapes Defined 


possible find the time peak, from Eq. 26. Proceed follows following 
Reich (9). 


Eq. can rewritten 


Thus, given and find from Eq. and then compute from Eq. 
28. The ratio m/G then found from Eq. using Fig. The hydrograph 
shape follows from Eq. 25. 
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APPLICATION 


For fully-gaged catchment, all the necessary parameters can esiimated 
described and the application straightforward. partial duration series 
floods used obtain the parameters, then the return period contours 
Fig. which also appears Table not necessarily give the return 
period years. get the desired return period years, T,, the selection 
the correct return period contours, use information from Table 
can obtained through the relationship given Eq. For convenience, Eqs. 


which represents the return period, years; represents the return 
period for the truncated distribution obtained from the partial duration series 
floods; and represents the mean number flood events per annum 
the truncated sample. find for the truncated distribution from the underlying 
total distribution, the following manipulation necessary. 

Let the return period for the total distribution the partial duration series 
estimated from the truncated distribution denoted The exceedance 


For the truncated distribution with truncation such that 
the area truncated then the cdf the truncated distribution given 


Use Eqs. and trial and error procedure will then allow the 
calculation the correct return period contour which must selected from 
Table order produce the family run hydrographs with the desired 
return period. 

approximate way obtain the desired return period contour comes from 
simply assuming that the number events below the truncation level 
proportional those counted above the truncation level. That means the number 
events expected per year over the full (untruncated) sample will 


HY6 
| 
approximation which will result the desired now simply 
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or, accurately 


that given and from the data, and the desired T,,, then the 
curve use. The following example should clarify the ideas. 
Example.—Assume that 50-yr return period family hydrographs 
desired for the Vaal River Standerton. From the streamflow record, the 
following information known: (1) Mean peaks, 6.162; (2) standard 
deviation peaks, 0.678; (3) mean volumes, 10.158; (4) 
standard deviation volumes, 0.907; (5) cross correlation coefficient, 
0.847; (6) area truncation, 0.447; (7) number events the 
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truncated sample 56; and (8) number years the record 54. 
Calculations.—The following calculations can made. 


The mean number events per year the truncated sample, 56/54 
1.037. 

Use Eq. gives the mean number events the underlying untruncated 
distribution: 1.037/(1 0.447) 1.875. 

The contour use Table therefore from Eq. 33: 1.875 

50.6 yr, which close enough the desired yr. 

Select pairs the standardized variables and from Table the 
contour (interpolate linearly for approximation between the contour and 


TABLE 3.—Parameters for Family Hydrographs with 50-Yr Return Period for Vaal 
River Standerton 


Hydrograph 
(selected meters per cubic meters 
line) second hours per second 


(4) 


the 100 contour) define suitable family run hydrographs. Since 
actual flood peak, exp cubic meters per second 35.31 
ft). For the actual flood volume: exp cubic meters per second 
hours hour 127,116 ft). 

The time from the peak the centroid the hydrograph given 
Eq. 27. 


now only remains solve for the time peak from and 
which Fig. useful. 

Each hydrograph now uniquely defined through Eq. 25. Otherwise, with 
the design hydrograph. 

family design hydrographs for the Vaal River Standerton with 
return period shown Fig. These hydrographs are from the 
parameters which appear Table Note that hydrographs and are from 
mirror images lines and Fig. 


HY6 
m/G, 
hours 
(1) (2) (3) (5) (6) (7) (8) 
1.325 2.284 55,800 2,231 4.20 1.8 
1.705 2.215 121,108 2,120 3.55 2.0 
2.030 2.030 162,199 1,878 4.11 2.8 
2.215 1.705 192,320 1,507 4.73 3.5 
2.284 1.325 204,740 1,164 5.33 4.5 
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Run hydrograph theory was developed and demonstrated for application 
fully-gaged catchments. The method allows for the synthesis the whole family 
design hydrographs for each desired return period. suggested that 
few members the family design hydrographs with the desired return period 
considered practical applications instead only the one member the 
family with the highest peak. This method should especially useful flood 
problems where the flood volume plays part. 

the development the theory, hydrographs were analyzed they occur 
the record; thus the sample hydrographs contains many composite hydro- 
graphs. This should result realistic design hydrographs for practical applications. 

the estimation the necessary parameters was demonstrated how 
straightforward can calculate the underlying probability distribution 
parameters from truncated sample. Truncation perhaps useful emphasize 
the effect the larger floods the distribution and ignore the smaller 
floods. 


The run hydrograph technique may also amenable towards extrapolation 
ungaged catchments. 
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The following symbols are used this paper: 


catchment area square kilometers 0.3861 mile); 
expected number events class interval 

probability distribution function for annual series floods; 
probability distribution function for partial duration series floods; 
probability density function; 

time from hydrograph peak its centroid, hours; 

number events; 

time from beginning flood hydrograph its peak, hours; 
number parameters estimated from sample; 

sample size; 

sample size for truncated sample; 

observed number events class interval 

discharge rate 35.31 cubic meters per second); 
peak discharge rate ft/sec; 

discharge rate ft/sec truncation level; 

coefficient association between variables; 

return period, years; 

return period (not necessarily years); 

return period approximation 

time; 

natural logarithm flood peaks; 

natural logarithm discharge truncation level; 

natural logarithm flood volumes; 

flood volume cubic meters 35.31 ft); 

flood volume cubic meters per second hours cubic meters 
per second hour 35.31 ft/sec hours); 

standardized natural logarithm for peaks, 
standardized natural logarithm for volumes 
function truncation 

estimate population mean flood; 

hydrograph characteristic which equals 

estimate population bivariate characteristic which 
gamma function; 

number degrees freedom; 

estimate truncation level pdf; 

estimate area pdf truncated; 
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mean number events per year truncated sample; 
mean number events per year untruncated sample; 
population mean; 

estimate population mean; 

moments round truncation level; 

population standard deviation; 

estimate population standard deviation; 

population cross-correlation coefficient; 

estimate population cross-correlation coefficient; 
standard normal density function; and 

chi-square. 
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DISCUSSION 


Note.—This paper part the Journal the Hydraulics Division, Proceedings 
the American Society Civil Engineers, Vol. 107, No. HY6, June, 1981. ISSN 
0044-796X /81 $01.00. 
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DISCUSSION 


The writers appreciate the interest shown the discussers. Some interesting 
points are addressed about the proposed model and dispersion theory general 
that need clarifying: the physical significance and time range values for the 
model parameters, the relative importance velocity distribution and storage 
mechanism, and the distinction between models that are time dependent and 
models with time dependent parameters. 

The parameter was defined measure the average number times 
particle goes into storage per unit time. The range the value 
shown Fig. over three log cycles, which large and justly criticized 
physical sense for this. However, this range for parameter evaluated 
for the cases reported Ref. 13. The range channel length and flow 
conditions for these cases shown Table The time range for then 
approximately the same the physical range for the data base. The ratio 
for each the cases (13) shown Table This ratio 
for particular case had range from 1.6 33.5 with average 7.2, 
opposed ratio over 1,000 for all data shown Fig. The variation 
exponentially decreasing function time that changes relatively 
quickly during the start dispersion and then undergoes only moderate change 
thereafter. 

The parameter measure the time particle moving downstream 
relative the total transit time. Holly and Tsai state that there apparent 
physical reason for this parameter small 0.75. Fig. shows the 
vertical velocity profile simulated the storage model. With this simplified 
velocity distribution, the parameter becomes lumped parameter that also 
takes into account the zone very low velocities along the bed and banks 
along with storage zones. Holly and Tsai also questioned the independence 
stage the value a,. Based limited data, was fairly constant 
for specific river reaches independent discharge and therefore stage. This 
may explained being indicative inherent roughness and storage 
zones which considered regards the total cross section may change 
proportionally with stage. The ratio for each the cases 
(13) shown Table 

The parameters vary systematically with time, however this does not change 
their physical significance. The initial phase dispersion occurs during the 
highest concentration tracer material, and there therefore high probability 
for transfer large mass dispersant into storage. The concentration 


“May, 1978, George Sabol and Carl Nordin, Jr. (Proc. Paper 13758). 

Assoc. Prof. Civ. Engrg., New Mexico State Univ., Box 3CE/Las Cruces, N.M. 
88003. 

U.S. Dept. the Interior Geological Survey, Box 15047, Water Resources 
Div., Stop 413, Denver Federal Center, Denver, Colo. 80225. 
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decreases with time and seems likely that the rate transfer into storage 
would also decrease with time. 

Prytch states that the writers used stochastic variables about which little 
information exists. The two parameters, and a,, are diffusion coefficients 
with real physical meaning; the fact that these coefficients are not the same 
those developed and investigated others should not important. While 
numerous studies have been undertaken document diffusion coefficients, 
the dispersion solution still eludes us. This suggests that the answer may 
found other approaches—divergent from the traditional. 

The relative importance storage mechanism dispersion difficult 
evaluate, and possibly dispersion studies are the most direct measure 
the importance the storage mechanism. The significance this storage may 
best illustrated the amount tracer material that recovered the 
most distant point downstream. The recovery ration, RR, defined the 
ratio mass tracer material measured site the mass tracer injected. 
The recovery ratio for each the most distant monitoring site for the cases 
(13) shown Table There are numerous difficulties and assumptions 
calculating recovery ratio and fact some recovery ratios are greater than 


TABLE 4.—Ranges Channel Length and Flow Conditions for Tracer Studies 
(13) Used Determine Model Parameters 


Variable 


(1) 


Minimum 


(2) 


Maximum 


(3) 


Length, miles 
Discharge, cubic feet per second 
Mean width, feet 
Mean depth, feet 
Mean velocity, feet per second 


one, however, the high incidence low recovery ratios indicates that storage 
mechanism may highly significant open-channel dispersion. 

The intent the proposed model was incorporate storage mechanism, 
not negate the importance the velocity distribution. The velocity distribution 
was simplified shown Fig. Prytch suggests that the numeric solution 
method Fischer (6) could used simulate these storage zones employing 
one more streamtubes zero velocity. that method the ratio zero 
velocity streamtubes total streamtubes would have determined, and 
method was suggested known for implementation this. 

has been generally argued for some time that the dispersion process 
open-channel flow procedes two phases. The initial phase, termed the 
convective advective period characterized incomplete mixing the 
tracer the cross section and skewed concentration distribution. The second 
phase, termed the diffusive, dispersive, Fickian period, obey the laws 
Fickian process given Eq. convenient and even physically appealing 
consider the dispersion process occur these two phases; however, 
reality doubtful that there distinct change between these, rather the 
change probably gradual. The dispersion process changes with time and either 
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TABLE 5.—Relative Mixing Time, Recovery Tracer Material, and Range Model 
Parameters for Tracer Studies (13) Used Determine Parameters 


(1) (2) (4) (5) 
40.0 0.32 10.8 1.46 
27.3 0.44 15.6 1.30 
59.4 0.50 4.1 1.35 
10.7 0.56 2.2 1.10 
3.5 1.22 3.6 1.09 
1.8 0.93 2.8 1.19 
3.8 1.25 2.1 1.09 
2.3 0.96 5.0 1.31 
3.2 0.44 2.6 1.21 
3.2 0.57 13.4 1.09 
2.8 0.93 4.0 
0.78 5.6 1.04 
9.9 0.93 6.3 1.08 
10.6 0.98 4.9 1.09 
4.1 0.51 2.8 1.09 
2.6 1.06 1.46 
48.6 0.65 4.4 1.22 
1.9 1.00 5.6 1.26 
18.4 0.61 1.6 1.09 
38.3 0.24 7.8 1.10 
67.6 1.40 19.2 1.21 
15.8 0.34 4.4 1.12 
55.4 0.78 6.9 1.18 
36.7 0.87 6.7 1.21 
17.4 0.70 3.6 1.20 
3.5 0.59 6.2 1.20 
3.3 0.77 4.1 1.30 
0.85 4.8 1.24 
1.5 0.62 3.7 1.15 
2.8 0.40 5.6 1.20 
5.0 1.00 33.5 1.42 
232.2 0.18 3.8 1.07 
110.6 0.16 2.7 1.15 
49.3 0.57 3.4 1.08 
0.9 0.64 5.2 
21.0 0.75 6.4 1.27 
23.8 0.52 7.8 1.26 
7.0 1.18 13.1 1.25 
1.7 1.03 9.4 1.22 
18.6 0.68 8.0 1.09 
4.1 0.71 2.1 1.12 
0.71 5.3 1.12 
4.5 0.86 3.5 1.32 
1.0 1.31 5.6 1.16 
1.9 0.85 28.4 1.14 
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TABLE 5.—Continued 


Maximum 
Minimum 

Mean 

Standard deviation 


Note: L/X total reach length mixing length; (RR) mass tracer recorded 
last site mass tracer injected, and and minimum values 
parameters for each case. 


the model the parameters should also change with time for adequate 
simulation the process. 

period time required the convective phase before the dispersant 
can considered completely mixed. This mixing length can approximated 


which the distance downstream from the source; the distance from 
the point maximum surface velocity the farthest bank; the hydraulic 
radius; and the shear velocity. The ratio reach length mixing length, 
for the cases given Table 

the cases reported (13) only two these took place entirely the 
convective period defined Eq. 32. the cases the minimum ratio 
reach length mixing length was 0.9, the maximum ratio was 232, and the 
average ratio was 18. The coefficient skew did decrease many these 
cases, but the distributions were still highly skewed the end the test reach 
even after having been the dispersive period for magnitude time larger 
than the convective period. Many the cases considered had been the 
dispersive phase for excessively long period and yet did not general behave 
according Fickian process. 

discussing the convective period, Prytch states that Eq. not appropriate 
during that time and that incorrect apply Eq. directly even after 
passing through the convective period; rather, the concentration distribution 
that exists the end the convective period must used the initial 
condition for the diffusive period. This distribution nor its location can 
accurately determined whereas the model presented Eq. valid 
throughout the complete dispersion process. 

Open-channel dispersion often considered one-dimensional Fickian process 
and solutions obtained using Eqs. and The application these 
equations however time dependent, that is, only valid after the initial 
convective period. Fickian theory, applied river dispersion, then itself 
time dependent. The proposed model can applied without consideration 


(1) (5) 
1.9 0.88 16.7 1.26 
1.3 1.13 13.3 1.13 
5.6 0.72 14.3 1.24 
232.2 1.40 33.5 1.46 
0.9 0.16 1.60 1.04 
18.7 0.743 7.198 1.186 
37.4 0.288 6.46 0.103 
Eq. 
R,u, 
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these convective dispersive periods; however, the parameters are time 
dependent. seems then that both procedures are fact time dependent; one 
virtue model applicability and the other virtue the parameters. 

Holly states that ‘‘the model should not expected applicable outside 
the ranges conditions and objectives for which was Table 
summarizes the range hydraulic conditions used forumulate the parameters 
the model. Applicability based this should quite high. 


Closure Chih Ted Yang’ and Charles Members, ASCE 


The writers appreciate the interests shown the paper. There dispute 
about the fact that solutions based the general theory minimum rate 
energy dissipation its simplified versions are good agreement with observed 
data from both laboratory flumes and natural rivers. However, seems that 
the discussers are not satisfied with some the assumptions and approaches 


used the writers their derivations. Most the disagreements between 
the discussers and the writers stem from different interpretations terms used 
the derivation and the concept used the paper which different from 
vectorial mechanics. general comparison between the vectorial approach and 
the variational approach will made clarify basic concepts. Specific questions 
raised each discusser will then replied separately. 

The science mechanics has been developed along two main lines, i.e., 
the vectorial and the variational approaches. The vectorial mechanics approach 
starts directly from Newton’s laws motion. Force and momentum are the 
basic concern the vectorial approach. Thus, the motion particle can 
uniquely determined the known forces acting every instant. The 
variational approach mechanics replaces the force the work the force 
the work function which frequently replaceable potential energy. 
the same time, the momentum Newton replaced kinetic energy. The 
variational approach scalar approach. this approach, the law mechanics 
stated the minimization the action integral. classical formulation, 
for nondissipative system, the action integral the integral kinetic energy 
minus potential energy. The mathematical tool used the vectorial approach 
the equation motion, while the calculus variation used solving 
minimization problems. Because the minimum value depends the constraints 

1979, Chih Ted Yang and Charles Song (Proc. Paper 14677). 

Engr., U.S. Dept. the Interior, Water and Power Kesources Service, Engrg. 


and Research Center, P.O. Box 25007, Building 67, Denver Federal Center, Denver, 
Colo. 80225. 


*Prof. Civ. Engrg., St. Anthony Falls Hydr. Lab., Univ. Minnesota, Mississippi 
River 3rd Ave., SE, Minneapolis, Minn. 
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which define the optimum path, the scalar approach that minimization 
with constraints instead solving differential equations with boundary conditions. 

Most the recent developments fluid mechanics and hydraulics are based 
the vectorial approach. This approach well suited solving problems 
which can described easily rectangular coordinate system, but becomes 
cumbersome for curvilinear coordinates. The process and resulting equations 
based the scalar approach remain valid for arbitrary choice coordinates. 
For systems with constraints, the scalar approach often simpler and more 
economical. Thus, the scalar minimization approach especially suitable 
for solving complicated problems natural rivers. more detailed discussion 
the two approaches and their merits given Lanczos (46). 

The general theory minimum rate energy dissipation states that when 
system equilibrium, its rate energy dissipation its minimum. 
This minimum value depends constraints applied the system. When the 
system not equilibrium condition, its rate energy dissipation not 
its minimum, but will decrease with respect time and reach minimum 
regain equilibrium. Depending the constraints and environment, there 
are different means which river can minimize its rate energy dissipation 
(48). The theoretical derivations made the writers the paper are limited 
the case that the only means available the system minimize the rate 
energy dissipation through the velocity distribution. The intention the 
writers’ derivations show that the variational minimization formulation 
the problem equivalent the use the equation motion for the case 
nonaccelerating open channel closed conduit flows. 

Silberman seems imply that because the developers turbulent theories 
the past did not use simple formula for the rate energy dissipation, 
the writers’ approach could not justified. The formulation Eq. based 
Boussinesq’s concept which makes possible define eddy viscosity for 
turbulent flows. Boussinesq’s concept assumption has been widely used 
turbulent theories. Whether this concept can directly used formulate the 
rate energy dissipation and the theory minimum rate energy dissipation 
depends whether the theoretical velocity profiles thus obtained agree with 
those observed turbulent flows. The writers’ verifications (41) indicate that 
the theoretical velocity profiles are good agreement with the measurements. 
The theory minimum unit stream power should applied one-dimensional 
open channel flows with constant width where sediment concentration not 
too high. Where the sediment concentration high, the power required transport 
sediment can longer ignored and treated constraint. this case, 
the rate energy dissipation for both water and sediment should considered. 
For most natural rivers with variable width, the theory minimum stream 
power, which also special case the general theory, should applied. 
Chang (43,44) applied the theory minimum stream power trial and error 
determine river pattern and regime channel geometry. His results agree very 
well with those observed the field. All these verifications (43,44) were made 
for turbulent flows natural rivers and regime channels. 

Davies considered laminar flow linear process; turbulent flow not. 
Actually, laminar linear nonlinear process depending whether 
the inertia term can neglected. the writers’ derivation, the eddy viscosity 
was not assumed constant, thus, the flow could nonlinear one. 
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Because represents the inertia effect turbulence and included the 
derivation, the inertia effect was not completely ignored. That is, the inertia 
effect the mean flow was ignored the derivation but not the inertia effect 
due turbulence. 

Parker raised several interesting points about the application the theory 
minimum stream power and minimum unit stream power. has indicated 
subsequent conversations with the writers that his discussion should not 
construed being blanket condemnation the approach power 
minimization. feels that worthwhile results have been obtained, but that 
the positive aspects some these might have been exaggerated. Parker’s 
discussion was addressed several specific points, and not the concept 
whole. The writers never stated that the theory minimum unit stream 
power equivalent the theory minimum stream power. However, both 
are special and simplified versions the general theory minimum rate 
energy dissipation. result using Eqs. and 30, Figs. and show 
that the measured values approach the theoretical values gradually with decreasing 
sediment concentration. When sediment concentration close zero when 
the bed incipient motion, the observed values are identical the theoretical 
values. The gradual transition agreement between observed and theoretical 
values function sediment concentration indication that the general 
theory could applied both rigid and movable boundary conditions and 
might bridge the gap between them. the determination velocity profiles 
open channel flows, the writers (41) treated the boundary rigid without 
sediment movement. For low medium sediment concentration, the case 
shown the original paper, sediment concentration was treated constraint 
the system. For the more general case, the total rate energy dissipation 
due water and sediment transportation may have considered (49). For 
all the aforementioned three categories verified the writers, observed values 
are good agreement with those determined from the theory minimum 
rate energy dissipation its simplified versions. Thus, the theory does have 
the potential bridge the gap between rigid and movable boundary fluid 
mechanics. 

The writers not believe that they have confused the minimization 
functional with the minimization functions. more detailed and rigorous 
treatment functional and function was made the writers recent paper 
(47). Theoretically speaking, the general theory minimum rate energy 
dissipation should applied wherever possible and the theory minimum 
stream power and the theory minimum unit stream power are special cases 
the general theory. These two simplified versions the general theory can 
applied under certain special conditions. Based comparisons with river 
and channel data, the theory minimum stream power can applied channels 
with variable width while the theory minimum unit stream power can 
applied channels with constant width. Among the sediment transport equations 
tested the writers, only Eq. can used determine the minimum unit 
stream power without using any resistance equation. Thus, Eq. can serve 
the dual purpose being transport well resistance equation. This 
was confirmed Parker shown Fig. Further studies the conditions 
under which simplified theories minimization can applied are needed. 

true that the use Eq. does not always assure the existence 
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minimum unit stream power for subcritical flows. shown Figs. and 
the agreements between observed and theoretical values deteriorate sediment 
concentration, which related Froude number, increases. For all the tests 
made the writers, minimum unit stream power can found supercritical 
flows. sediment concentration Froude number increased certain value, 
the treatment sediment concentration merely constraint the minimization 
process not sufficient. this case, the total rate energy dissipation for 
water and sediment has minimized (49). 

Chen raised some basic questions the application the variational principle 
turbulent flows. detailed and more rigorous treatment this subject was 
published the writers (47). The writers agree with Chen that the term 
varied flow’’ used the writers could cause some confusion and 
misunderstanding. The ‘‘gradually varied flow’’ used the paper means nonac- 
celerating flow the inertia terms the equation motion are small and 
can ignored (47). The writers not believe that the dissipation function 
for turbulent flow was incorrectly expressed Eqs. and 7b. There are 
least two ways demonstrate this point. The simplest way start with 
the Reynold’s equation, Eq. and the Boussinesq’s assumption, Eq. Substitu- 
tion Eq. into Eq. leads equation the form identical that 
the original Navier-Stokes equation based mean values with replacing 
Eq. simplified version that equation after dropping the acceleration 


terms. The dissipation function based the instantaneous values known 
be: 


au, 


Eq. the dissipation function the flow represented the Navier-Stokes 
equation and the reduced equation motion for turbulent flow obtainable 
from the Navier-Stokes equation substituting the instantaneous values with 
the mean values and the same substitution, Eq. transformed 
into Eq. 7a. That is, Eq. the dissipation function the Navier-Stokes 
equation, then Eq. must the dissipation function the flow represented 
Eq. There should doubt that one can follow, step step, the 
same procedure required for the derivation Eq. arrive Eq. starting 
from Eq. 

The second approach base the argument the energy balance equation 
shown Hinze (45). Chen has correctly pointed out, the direct viscous 
dissipation kinetic energy into heat energy from the mean velocity field and 
the turbulence field are, respectively, given Eqs. and 36. However, 
Hinze (45) and other investigators have pointed out, the process energy 
transfer within unit volume fluid is: mean flow large scale turbulence 
small scale turbulence heat energy. The rate turbulence production 
the expense mean flow energy per unit volume is: 


Boussinesq’s assumption, Eq. reduces the second part Eq. 
pointed out Chen. That is, Eq. equivalent to: 


(44) 
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Now all that remains shown that under equilibrium condition. 
Clearly, this not strictly true because the presence other nondissipative 
work done terms. This strictly true steady uniform turbulent flow because 
this case there net work done between neighboring volumes fluid 
and the turbulent energy each unit volume remains constant. That is, the 
turbulence production strictly balanced with the turbulence dissipation every- 
where this case. For the nonaccelerating flow case treated the writers’ 
paper, the relationship only approximately true. When the system whole 
equilibrium condition, the only source error the differential turbulent 
energy flux and the differential work done between the upstream and downsiream 
ends. For very long and nearly nonaccelerating river, these differences are 
likely much less than the total energy loss. The work done the solid 
boundary well the free surface boundary zero because either the velocity 
the stress zero there. The writers feel that Eq. reasonably good 
representation the total energy dissipation rate for the intended type flows. 

appears that Chen might have misinterpreted the meaning One must 
variation and not the velocity. Therefore, the velocity distributions both 
ends study reach are given, then must zero there. Because the free 
surface stream surface where the mean velocity must zero 
along the stream surface. the velocity distributions both ends are not 
given, the hydraulic heads the upstream and downstream ends study 
reach the difference between them must given before the problem can 
solved (47). The writers agree with Chen that the rigid boundary gradually 
varied flow problems given discharge with given roughness and specified 
boundary conditions can solved without resorting minimization approach. 
The intent the writers’ derivations was show that the variational 
minimization approach can provide alternative, and some cases more 
desirable and flexible approach solving open channel hydraulic problems. 
The basic concept behind the minimization approach not limited seeking 
the best estimated value roughness coefficient. 
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INTERFACIAL STABILITY CHANNEL 


The writer appreciates the interest shown the paper and wishes thank 
the discussers, Jensen and McCutcheon, for their comments and suggestions. 
With regard these discussions, the writer would like note the following. 

First, the suggestion Jensen that interfacial stability large Reynolds 
numbers should function and would appear valid. Jensen 
essentially claiming that interfacial stability function the bottom boundary 
friction factor which not violation the writer’s basic hypothesis. 

Second, McCutcheon correctly noted that the aspect ratio (width-to-depth) 
the flume used the writer was very small. McCutcheon (19) has demonstrated 
that low aspect ratio can have significant effect stratified flows. 
This discusser’s comments would indicate that the corrective procedures devel- 
oped Vanoni and Brooks (17) are not entirely applicable stratified flows. 

Third, McCutcheon noted that the writer’s hypothesis showed improvement 
the work Polk al. (12) with regard upstream warm water wedges 
power plants. Although the writer noted and commented the work 
Polk, al. (12), the specific subject this paper was the stability superposed 
turbulent layers with zero velocity difference the interface. this case, 
there single mechanism turbulence generation; the bottom boundary. 
the case upstream warm water wedges, there are two mechanisms 
turbulence generation; the bottom boundary and interfacial shear. Because 
this difference turbulence generation, these two cases should considered 
separately. 

Fourth, the writer would like thank McCutcheon for publishing additional 
data which were not available when the writer’s work was published. addition, 
the error analysis done this discusser interesting. 

conclusion, the writer would note the recent work Abraham, al. 
(24) the areas interfacial shear and stability. 


24. Abraham, G., Karelse, M., and Van Os, G., the Magnitude Interfacial 
Shear Subcritical Stratified Flows Relation with Interfacial Journal 
Hydraulic Research, Vol. 17, No. 1979, pp. 273-288. 


“August, 1979, Richard French (Proc. Paper 14768). 
Assoc. Research Prof., Desert Research Inst., Water Resources Center, Univ. Nevada 
System, 4582 Maryland Parkway, Las Vegas, Nev. 89109. 
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UNCONDITIONAL STABILITY CONVECTION 


Closure Victor Miguel Ponce,” Yung Hai Members, ASCE, 
and Daryl ASCE 


The writers wish thank Stelling, Arunachalam and Raman, and Barnett 
and McBride for their interest this paper. 

Stelling has addressed the question terminology explicit versus implicit 
schemes the one-dimensional convection equation. addition, has pointed 
out some the complications which may arise calculations involving tidal 
areas. Arunachalam and Raman’s discussion centers the similarity (or perhaps, 
dissimilarity) the writers’ equations with the well-known Crank-Nicolson 
scheme. Barnett and McBride have correctly pointed out some the difficulties 
encountered when attempting model functions containing sharp gradients using 
neutrally stable schemes. 

The semantic question regarding explicit and implicit schemes one that, 
the writers’ opinion, needs put rest. has long been established 
(4) that explicit means point-by-point calculation; implicit, the other hand, 
means line-by-line calculation. this terminology, the schemes treated the 
writers are all explicit. particular application the boundary conditions may 
such that the solution necessitates the inversion matrix, which case 
the associated scheme would implicit. However, the basic difference between 
explicit and implicit schemes remains intact. 

Arunachalam and Raman through lengthy and rather pointless mathematical 
exercise attempt show that the writers’ equations bear certain similarity 
with the Crank-Nicolson scheme. fact, the writers treat explicit numerical 
schemes the pure convection equation; Crank-Nicolson’s method well- 
known implicit scheme the diffusion equation. The only commonality between 
the writers’ approach and that Crank-Nicolson seems the feature 
space and time centering. 

Barnett and McBride’s comments are well taken. The writers’ aim was 
show the numerical properties class explicit schemes the one-dimensional 
convection equation. doing so, they may have inadvertedly misled the reader 
into thinking that schemes with vanishing numerical diffusion coefficient are 
preferred all cases. fact, Barnett and McBride point out, there 
something say favor the introduction some numerical diffusion, 
long does not considerably impair accuracy. 

The writers did not intend use the schemes presented connection with 

1979, Victor Miguel Ponce, Yung Hai Chen, and Daryl Simons 
(Proc. Paper 14807). 

Prof. Civ. Engrg., Coll. Engrg., San Diego State Univ., San Diego, 
Calif. 92182. 

Assoc. Prof. Civ. Engrg., Colorado State Univ., Fort Collins, Colo. 80521. 

Dean Engrg. and Prof. Civ. Engrg., Colorado State Univ., Fort Collins, 
Colo. 80521. 
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the computation These waves contain, definition, very 
sudden changes properties, and numerical method can fully account for 
these features. Barnett and McBride’s choice minimal phase error being 
more desirable than amplitude preservation best compromise, given the 
limitations inherent the discrete computation. 


The discussion Oner Yucel requests additional clarification the nature 
approximations applied hydraulic graduates and flowrates which allows 
linearization the nonlinear network equations. There are, fact, approx- 
imations involved because essentially exact solution the nonlinear network 
equations are obtained. Linearization the nonlinear equations utilized 
iterative procedure for solving these equations, but carrying out several trial 
results the true solution the nonlinear set equations. 

Certainly the ultimate goal cost optimization slurry transport systems 
desirable and the technique described the paper for handling the hydraulic 
calculations essential component this type study. 


Discussion Spyridon 


The writer agrees with the main thesis the paper, stated the following 
sentence: effects dead zones well channel irregularity and 
nonuniformity can regarded mere prolongation the time required 
reach the Taylor diffusive Using this same postulate, the writer 
(1,30) arrived different route variation with which equivalent 
that the authors. Putting and using the approximate transformation 


1980, Don Wood (Proc. Paper 15132). 

Civ. Engrg., Coll. Engrg., Univ. Kentucky, Lexington, Ky. 40506. 

1980, Henry Liu and Alexander Cheng (Proc. Paper 15454). 

Scientist, Hydr. Div., National Water Research Inst., Canada Centre for 
Inland Waters, P.O. Box 5050, Burlington, Ontario L7R 4A6, Canada. 
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(also used the authors implicitly), Eq. reduces Eq. Ref. 
30; note that the writer’s (30) coefficient confused with the 
function the authors (therefore quotation marks are used this discus- 
sion)—is equal LU. may interest compare next the experimental 
values corresponding model parameters. 

Considering Eq. and putting 0.7 L/U gives: 


which coincides with Eq. Ref. 0.5 set equal the writer’s 
coefficient a,. Using recommended the authors, gives 
2.5; this about twice the average value implied Fig. Ref. but 
still within the range scatter the data points. Rearranging Eq. with 


= 
RU, 
which coincides with Eq. Ref. set equal the writer’s coefficient 
The authors found vary from 0.2 1.3, having average value 
about 0.6. These are close corresponding values implied Fig. 


Ref. 30, though seems scatter more than Considering the writer’s coefficient 
Eqs. and may used simultaneously give 


U 2 


which agrees closely with average values implied Fig. Ref. 
30. summary, may stated that the experimental values model parameters 
determined the authors and the writer, respectively, are agreement. 
general, the writer’s ranges scatter appear wider than those the authors; 
this probably due the fact that the authors have allowed larger scatter 
evaluating individual test results can seen comparing Figs. and 


Concerning details the model proposed the authors, the writer suggests 
the following: 


Introduction the Pearson distribution (Eq. 20) predict 
versus variations appears arbitrary, especially view the fact that the 
model does provide pertinent equation (Eq. 10). The advantage cited the 
authors for Eq. 20, i.e., that reproduces the skew better than does Eq. 10, 
considered marginal from the practical point view. practice, the errors 
committed due uncertainties predicting and far outweigh discrepan- 
cies due inaccurate predictions skew. Therefore, use Eq. considered 
unnecessary refinement. 
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direct substitution and some algebra, can shown that Eq. 
not solution Eq. 11; the correct version Eq. 11, such that satisfied 
Eq. 10, can shown be: 


Note that, Eq. 0.5 which consistent with the 
definition the diffusion dispersion coefficient; and that does 
not coincide with unless constant. According Eq. 8a, function 
approaching the constant value only when Using the 
definition and performing the differentiation indicated Eq. gives: 


which shows that the dispersion coefficient, increases with time 
the asymptotic value For practical purposes, the dispersion coefficient 
can considered equal when The correct version Eq. 


ac’ 


noted, however, that these modifications not affect later results because 
the analysis has been based and which are correct. 


30. Beltaos, S., Dispersion Journal the Hydraulics Division, 
ASCE, Vol. 106, No. Proc. Paper 15118, Jan., 1980, pp. 151-172. 


The writer welcomes the new experimental data given this paper, accurate 
observations transverse mixing real streams are indeed few. Nevertheless 
the writer concerned that the author’s computed results are strongly 
dependent the factor and would like more certain that the reported 
values the transverse mixing coefficient are accurate. earlier work 

“October, 1980, Spyridon Beltaos (Proc. Paper 15731). 


Civ. Engrg., Hydr. and Coastal Engrg., 412 O’Brien Hall, Univ. California, 
Berkeley, Calif. 94720. 
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(16) the author has also used the stream tube method ascertain the value 
but did not require the assumption that constant. Since 
proportional which far from constant. addition, should 
noted that Eq. follows from Eq. only constant, since otherwise 
the averaging process yields additional term resulting from the correlation 
with U,. The definition (given incorrectly Eq. and the derivation 
Eq. depend turn Eq. Eq. being the crucial equation for the 
determination from the observations. Without access the complete 
set data the writer unable evaluate the inaccuracy introduced the 
several averaging processes. The author notes briefly that one instance 
checked his value using numerical algorithm. Any further evidence 
that could offered the closure establish the validity the reported 
values would most helpful. Use the numerical algorithm described 
the writer Ref. 16, e.g., would avoid the questions accuracy just raised. 

The writer hopes that maps the test reaches for tests 3-4 and 5-6 can 
included the closure, suspects few readers will have access 
Ref. and thinks maps the study reaches essential evaluate results. 

Finally, the writer wishes suggest explanation for the result that 
K(=E,/RV,) increases due ice cover. wide stream and 
the nondimensional representation commonly used since Elder’s 
(1959) work. wide ice-covered stream, however, H/2. The question 
then arises whether the correct length scale, i.e., whether more 
E./HV, K/2. the writer’s opinion the choice clear. Taylor’s 
theory diffusion continuous movements (21) shows that the magnitude 
turbulent diffusion coefficient depends the product length scale 
and velocity scale (see Ref. 20, 71). turbulent channel flow the velocity 
scale generally accepted the shear velocity. The length scale generally 
taken the distance between the boundaries, which means between the free 
surface and the channel bottom the case free surface flow, between 
the upper ice boundary and the channel bottom the case ice-covered 
stream. The author’s use implies that the effect the ice cover 
halve the length scale the turbulence. This seems unlikely, since the turbulent 
motions can expected persist over the whole distance between the 
boundaries, and the writer suggests that the logical choice. 

Table gives information computed from the author’s Table but the 
nondimensional result given terms E./HV,. The shear velocity, 
was computed from the author’s table assuming, for lack 
more exact information, that for open water and for ice-covered 
conditions. 

The range nearly within that previous experiments 


reviewed the writer (Ref. 20, 112). two the three pairs, 
smaller for the ice-covered test, but consistent relationship shown. 
The writer concludes that the author’s data support the previously established 
approximate relationship, and demonstrate that ice cover has great effect. 
the same time the author’s data demonstrate once again our inadequate 
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TABLE 3.—Recomputed Dimensionless Dispersion 


Discharge, 
Riversite cubic 
and channel Flow meters Width, Depth, 
description condition per second meters meters 
(2) (3) (4) (5) (6) 

Athabasca River Ice-covered 

below Fort Mc- 

Murray; straight 

with occasional 

bars, islands 
Same Test Open water 
Beaver River near Open water 

Cold Lake; regular 

meanders, point 

bars and large 

dunes 
Same Test Ice-covered 
Athabasca River Open water 

below Athabasca; 

irregular meanders 

with occasional 

bars, islands 
Same Test Ice-covered 


Notes: 3.28 ft. Col. 1-6, and are taken directly from Table 


understanding the nature transverse mixing, exemplified the range 
the dimensionless coefficient over factor 


19. Elder, W., dispersion marked fluid particles turbulent shear 
Fluid Mech. 544-460 (1959). 

20. Fischer, B., List, Koh, Y., Imberger, J., and Brooks, H., Mixing 
Inland and Coastal Waters, Academic Press, New York, N.Y., 1979. 

21. Taylor, Continuous Proceedings, London Mathe- 
matical Society, Series A., Vol. 20, 1921, pp. 196-211. 
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for Tests 


Transverse 
mixing 
Shear coefficient, 
Hydraulic Velocity, velocity, square 
radius, meters meters meters 
meters per second per per second 


(7) (8) (10) (11) 


Discussion Peter and Lam 


The authors have presented method for predicting the stage discharge for 
gravel-bed streams. Their contribution has been greatly enhanced the use 
relatively rare and comprehensive data from actual river reaches. Fundamental 
the authors’ analysis was the determination the bar resistance and 
this aspect their paper that the writers wish comment. 

“October, 1980, Gary Parker and Allan Peterson (Proc. Paper 15733). 


Engr., Environmental Hydr. Section, Hydr. Div., National Water Research 
Inst., Canada Centre for Inland Waters, P.O. Box 5050, Burlington, Ontario L7R 4A6, 
Canada. 

Environmental Hydr. Section, Hydr. Div., National Water Research Inst., Canada 
Centre for Inland Waters, P.O. Box 5050, Burlington, Ontario L7R 4A6, Canada. 


795 
(12) (13) 
0.95 0.49 1.44 0.0366 0.041 1.18 0.59 
2.2 0.95 1.44 0.0557 0.093 0.76 0.76 
0.96 0.50 2.1 0.0444 0.043 1.01 1.01 
0.305 0.28 2.1 0.025 0.020 2.62 1.31 
2.05 0.86 3.1 0.0789 0.067 0.41 0.41 
0.48 0.40 3.1 0.0382 0.010 0.55 0.27 
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One the two fundamental assumptions made the authors that 
unique function the dimensionless shear stress due the granular 
bed material. The values obtained were plotted against Fig. 
arrive relationship between and The writers wish suggest 
simpler variable which may correlate just well with C,. 

The submerged specific gravity will neglected was taken constant. 
Following Keulegan’s relation for grain resistance the authors obtained Eq. 
which gives function and the Froude number Therefore, 
the authors are effect assuming that function variable which 


with given Eq. 10. 

The fact that should depend and can obtained from 
dimensional analysis for and then neglecting the effects Reynolds number 
and the specific gravity the bed material. However, has been shown 
Engel and Lau (35) that, for Froude numbers 0.5, the total resistance 
for flow over dunes was independent the Froude number. the grain resistance 
coefficient also independent the Froude number (Eq. 2), one would expect 
that, the present case, the bar resistance coefficient will also independent 
function only i.e.: 


The authors’ Fig. for Old Man River shows that even though the effect 
not too large, has large effect the total and 
therefore C,. Therefore, may better try and correlate with 
which more basic variable opposed the more complicated This 
would also more directly reflect the authors’ premise that becomes more 
important successively lower stages the flow. 

second point which the writers wish comment the authors’ assumption 
that, flood stages, negligible. The authors sought confirm their 
assumption comparison and S,, their Fig. and observed that 
the which they calculated was approximately equal the observed 
flood stages. These results are contrast with observations made artificial 
dunes the writers from laboratory measurements (35). Ref. the 
writers’ experimental data, well those Vadagodakumbura for values 
1,000 and dune steepness (i.e., dune length) down 
0.02, indicate clearly that the total friction factor larger than that computed 
from the authors’ Eq. may that, because errors determining the 
slope and grain-size the field, the difference between and S,, could 
not detected the authors. seems only logical that even high flows 
there must some form resistance due the shape the bars. Therefore, 
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depending upon the shape the bars, may small but may not always 
negligible. 


35. Engel, and Lau, L., Factor for Two-Dimensional Dune Roughness,”’ 
Journal Hydraulic Research, Delft, The Netherlands, Vol. 18, No. 1980, pp. 
213-225. 
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